# **SpringerBriefs in Computer Science**

**David Ryckelynck · Fabien Casenave · Nissrine Akkari**

# **Manifold Learning** Model Reduction in Engineering

**SpringerBriefs in Computer Science**

SpringerBriefs present concise summaries of cutting-edge research and practical applications across a wide spectrum of fields. Featuring compact volumes of 50 to 125 pages, the series covers a range of content from professional to academic.

Typical topics might include:


Briefs allow authors to present their ideas and readers to absorb them with minimal time investment. Briefs will be published as part of Springer's eBook collection, with millions of users worldwide. In addition, Briefs will be available for individual print and electronic purchase. Briefs are characterized by fast, global electronic dissemination, standard publishing contracts, easy-to-use manuscript preparation and formatting guidelines, and expedited production schedules. We aim for publication 8–12 weeks after acceptance. Both solicited and unsolicited manuscripts are considered for publication in this series.

\*\*Indexing: This series is indexed in Scopus, Ei-Compendex, and zbMATH \*\*

David Ryckelynck · Fabien Casenave · Nissrine Akkari

# Manifold Learning

Model Reduction in Engineering

David Ryckelynck Centre des Matériaux Mines Paris—PSL Évry, France

Nissrine Akkari Department of Digital Sciences and Technologies Safran Tech Châteaufort, France

Fabien Casenave Department of Digital Sciences and Technologies Safran Tech Châteaufort, France

ISSN 2191-5768 ISSN 2191-5776 (electronic) SpringerBriefs in Computer Science ISBN 978-3-031-52766-1 ISBN 978-3-031-52764-7 (eBook) https://doi.org/10.1007/978-3-031-52764-7

This research was partially funded by the French Fonds Unique Interministériel (MOR\_DICUS).

© The Editor(s) (if applicable) and The Author(s) 2024. This book is an Open Access publication.

**Open Access** This book is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.

The images or other third party material in this book are included in the book's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the book's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

The use of general descriptive names, registered names, trademarks, service marks, etc. in this publication does not imply, even in the absence of a specific statement, that such names are exempt from the relevant protective laws and regulations and therefore free for general use.

The publisher, the authors, and the editors are safe to assume that the advice and information in this book are believed to be true and accurate at the date of publication. Neither the publisher nor the authors or the editors give a warranty, expressed or implied, with respect to the material contained herein or for any errors or omissions that may have been made. The publisher remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

This Springer imprint is published by the registered company Springer Nature Switzerland AG The registered company address is: Gewerbestrasse 11, 6330 Cham, Switzerland

Paper in this product is recyclable.

### **Foreword**

Safran Tech is the corporate research center of Safran Group, a French multinational company that designs, develops and manufactures aircraft and rocket engines, as well as various aerospace and defense-related equipment. This book is the result of the collaboration of two of my research engineers at SafranTech, Fabien Casenave and Nissrine Akkari, with David Ryckelynck, professor of Applied Mathematics at MINES Paris PSL, a French engineering school. This collaboration dates back to the creation of Safran Tech in 2015 and has involved Ph.D. students, namely Thomas Daniel and Hamza Boukraichi. Thomas Daniel won the 2022 AMIES Ph.D. award (AMIES: Agency for Mathematics in collaboration with companies and the society) for his work on dictionary-based ROM-nets, which are introduced in this book.

Nowadays, the engineering and design processes involved when creating a new mechanical part, component or product in industrial companies like Safran rely on some level of numerical simulation. In some cases, complex and high-dimensional models are required to be computed, and often in a so-called "many queries" context: when searching for an optimal design, quantifying uncertainties or exploiting a "realtime" simulator, the physics problem is solved a large number of time at different configurations. Hence, we need fast surrogates, with a control of the approximation error. Among the available families of methods, this book deals with *projectionbased reduced-order models*, which consists in solving the same physics equations as the reference high-dimensional models, but on a reduced dimensional vector space. Resorting to the known physics equations provides explainability, error estimates in some cases, and good performance with sparse learning databases. Classical and modern (e.g., deep learning) machine learning technologies are used to assist *projection-based reduced-order models*, either in annex tasks like clustering or classification or by providing non-linear dimensionality reduction procedures. Challenges of these methods include geometrical variabilities, predictive uncertainty quantification, and efficiency (i.e., speed-up factor, also taking the learning stage into account).

This book is addressed to graduate students and researchers who look for an introduction to *projection-based reduced-order models*, an understanding of the main features and challenges of these methods, as well as its state of the art application to real-life engineering design problems and its extensions.

Paris, France September 2023 Christian Rey

**Acknowledgements** We acknowledge the contribution from the many co-authors of our publications, in particular Thomas Daniel and Christian Rey for their help in Chap. 5.

**Christian Rey** is head of the *Mathematics and Algorithms for Design by Simulation* research unit at SafranTech and professor of Computational Mechanics. Christian is a long-time collaborator of the authors of this book.

## **Contents**



### **Acronyms**



## **Chapter 1 Structured Data and Knowledge in Model-Based Engineering**

### **1.1 Nomenclature**

In the book, the following notation is used: bold face symbols denote vectors (lowercase letters) . **a**, . **b** ... or matrices (uppercase letters) . **A**, . **B**... Data are denoted by .**X** for input data and labels or output data are denoted by . **Y**. The spatial gradient operator is denoted by. ∇, and the divergence is expressed by.∇·•. The dependence on spatial coordinates. *ξ* is most of the time omitted for simplicity of notation. Modeling parameters are denoted by . *μ* or .*μ*□ when they have a tensor shape (multiple indices). Weights in neural networks are denoted by .**W**. In mechanical models, the Cauchy stress is denoted by.*σc*. Solutions of physics-based models are denoted by. ~*<sup>u</sup>* for partial differential equations (PDEs) or ordinary differential equations (ODEs). High-fidelity solutions in finite dimensional space are denoted by . *u*. Solutions in reduced dimensional spaces are denoted by. *u*. PDEs are set up over a domain denoted by. Ω. Solution spaces are the following:


The finite element shape functions are denoted by . ϕ. The coordinates of a finite element solution are denoted by the vector. **u**. The reduced vector space.V*<sup>n</sup>* is spanned by empirical modes denoted by .(ψ*<sup>k</sup>* )*<sup>k</sup>*=1,...,*<sup>n</sup>*. The matrix form of a reduced basis is denoted by .**<sup>V</sup>** <sup>∈</sup> <sup>R</sup>*<sup>N</sup>*×*<sup>n</sup>*. The vector of reduced coordinates is denoted by .*<sup>γ</sup>* <sup>∈</sup> <sup>R</sup>*<sup>n</sup>*; . γ*<sup>k</sup>* ,.*k* ∈ [1,..., *n*] is the.*k*th component of. *γ* . When needed, variables depending on parameter. *μ* are denoted using an exponent.•*<sup>μ</sup>* or an index.•*μ*.

### **1.2 Model-Based Engineering**

Today, a large set of engineering tasks is supported by mathematical models related to various scientific disciplines. This set of tasks is called **Model-Based Engineering** (MBE). In this book, we restrict our attention to geometrical or morphological models, thermal models, mechanical models and statistical models, at various scales. The related mathematical equations can be algebraic equations, ordinary differential equations, partial differential equations (PDEs), or every possible combination of these types of equations in a coupled system of equations.

Various numerical methods have been developed in applied mathematics, to find approximate solutions of these different types of equations. Therefore, mathematical properties are available related to the convergence of the numerical solutions to the exact solutions having no numerical approximation. The engineers have a strong confidence in the numerical predictions obtained by these models, provided that they are well used in their domain of validity. Hence, in practice, numerical models and numerical predictions are ubiquitous in MBE.

In industrial sectors, such as automotive industry, aeronautical industry, energy, shipbuilding, electronics manufacturing, etc., **computer platforms for MBE** have been developed for assemblies of components, to ease the communication of data and models, between engineering tasks. Quality and efficiency in engineering tasks depend directly on the **continuity of data** on computer platforms. As explained in the ATLAS program of AFNet, <sup>1</sup> such continuity of data requires the development of standards for data exchanges between stakeholders in a project, a production line, a department, etc. Assemblies of components that are equipped with numerical models are termed complex systems in this book.

A **complex system** is an assembly of components that have a numerical representation related to coupled physics-based models, for its design phase, its manufacturing phase or its exploitation phase.

In computer platforms for MBE, complex systems have an idealized digital representation, termed **digital mock-up**, where model couplings can be implemented in a numerical model dedicated to an engineering task. For instance, in continuum mechanics, the mechanical model of a component is weakly coupled to the geometrical model of this component. Many multi-physics couplings in coupled problems, can be implemented on computer platforms for MBE. Complex numerical simulations in MBE are supported by computers that range form laptops to computer clusters for high-performance computing (HPC). In most cases, geometrical models are created by using computer-aided design (CAD) in order to obtain a common geometric model for both design tasks and manufacturing tasks, of each part and assemblies in a complex system.

Data flux in computer platforms for MBE are related to design tasks, manufacturing tasks, non destructive testing (NDT) tasks, exploitation or monitoring tasks, inspection tasks, repairing tasks, recycling tasks. Figure 1.1 is a summary of these data fluxes.

<sup>1</sup> https://www.afnet.fr/en/.

**Fig. 1.1** Schematic view of data continuity in an computer platform for Model-Based Engineering (MBE), enabling engineering tasks related to CAD, PDE solution

The value of data, beyond the tasks at hand, should not be ignored by engineers. In particular, the value of data must be considered in light of advances in artificial intelligence (AI) and machine learning methods.

### **1.3 Motivation**

Our main paradigm is: the more we know about underlying mathematical properties or physical properties of data, the less we need data for machine learning tasks. In this book, we aim to learn manifold from data and conventional knowledge in engineering. This paradigm has been recently employed in Physics Informed Neural Networks [ 6, 7] and knowledge transfer [ 4]. Similarly to transfer learning, we are "motivated by the fact that people can intelligently apply knowledge learned previously to solve new problems faster or with better solutions" [ 5].

Much data in computer platforms for MBE are structured data, in the sense that their format is predefined and documented (this notion is different from structured and unstructured meshes, referring respectively to constant rectilinear and heterogeneous meshes). These structured data are supported by a graph or they have a tensor format, like digital images. Since deep learning has made huge progresses in computer vision, image transformation and for graph neural networks, there is an opportunity to develop deep learning in MBE. The development of sensors and image acquisitions in non destructive testing of complex systems makes this opportunity more and more relevant. The connection of computer platforms to HPC facilities leads to the two following issues:


In this book, we propose numerical methods and algorithms to assist engineers in various conventional tasks by using machine learning. This gives value to data that flow in computer platforms for MBE. Our goal is to augment engineers capabilities by using artificial intelligence, without ignoring the knowledge already available for engineering in all scientific disciplines. We think that the collaboration between engineers and data scientists is going to grow up, but engineers will not be replaced by them. At the engineering school MINES Paris PSL, professors David Ryckelynck and Elie Hachem created in 2017 master courses to ease such collaborations between engineers and data scientists. The title of this set of courses is Ingénierie Digitale des Systèmes complexes (Digital engineering of complex systems). The creation of such a course would not have been possible without the help of SAFRAN, both for the illustrations of classical AI and for the development of an engineering augmented by a hybrid AI, which preserves scientific knowledge.

Acceleration of engineering tasks is the main objective here, but we are also able to develop engineering tasks that were unaffordable without machine learning. This is especially true when considering uncertainty quantification (UQ) or image-based digital twining of complex systems.

Following the definition in [ 2], a digital twin involves similar numerical models to mock-up but it is related to an existing instance of a complex system or an instance of one of its component. A digital twin also include the related observational data. Hence, in essence, a digital twin has **multimodal data**, such as observational data and simulated data. The acceleration of image-based modeling for digital twining of manufactured parts is going to foster the development of functional assessment, in which defects harmfulness are evaluated with respect to physics-based of mechanicalbased specifications.

The knowledge accumulated via digital twining or via sensors only, enable the development of uncertainty quantification. Sensors and numerical mock-up being more and more spatially resolved and accurate, the acceleration of numerical predictions is becoming crucial for practical application in engineering. In Sect. 5, mechanical predictions are accelerated by using machine learning in order to account for stochastic variation of temperature fields when computing the lifetime of turbine blades.

In this book we develop the following machine learning tasks, via manifold learning:


Data visualization is out of the scope of this book, although dedicated manifold learning techniques have been proposed in the literature. In this book, manifold learning is thought of as an attempt to generalize linear frameworks for numerical methods in engineering.

When accumulating data in computer platform for MBE, data do not fill the entire ambient space used for the encoding of these data. There exists an empirical space, or a **latent space**, occupied by data whose dimension is usually much smaller than the ambient space. Dimensionality reduction and manifold learning aim to find a numerical approximation of this latent space. As models are ubiquitous, we have to consider large sets of synthetic data, or more precisely, sets of simulated data for engineering tasks. These simulated data are solutions of numerical equations. They belong to a **solution space**, which is the ambient space for simulated data. A solution space is directly related to the general numerical scheme used to compute instances of numerical predictions. Given an engineering task, an ambient space is a common space for all data instances. It does not depend on any variable parameter. The dimension of the ambient space is set up so that each instance of data is a point inside the ambient space.

Two kinds of latent spaces can be defined when considering simulated data:


Projection-based model order reduction is detailed in Chap. 2. Our recommendation for selecting a surrogate model rather than a projection-based reduced order model is the following. If you consider a mature engineering task so that:


then you should use a surrogate modeling rather a projection-based reduced order modeling. But, when modeling a complex system by using PDEs, the selection of simulation outcomes is rarely trivial. Generally speaking, projection-based reduced order models are slower than surrogate models, but one can adapt the outcomes of reduced predictions and their validity domain is not restricted to the parameter domain used to learn the reduced order model, both in term of extent and in term of dimension.

In practice, data accumulation has to be supplemented by data compression schemes but also data pruning scheme. As explained in [ 3], manifold learning for projection-based model order reduction is a convenient way for data pruning that preserves modeling capabilities. Pruning algorithms facilitate data augmentation for high-dimensional data, by limiting the memory requirements for augmented data [ 1].

So, data in MBE are simulated data and observational data, where simulated data are synthetic data forecast by physics-based models. Unfortunately, observational data on complex systems are extremely rare in engineering. For instance, aircraft landing gear tests, automotive crash tests, etc. are mostly done for conformity assessment processes or for certification processes. For example, in the aeronautics industry, aircraft manufacturers must demonstrate compliance of new products with regulatory requirements. This compliance demonstration is done by analysis during ground testing (such as tests on the structure to withstand bird strikes, fatigue tests and tests in simulators) but also by means of tests during flight. Usually, manufacturers have many more instances of observational data, at the components scale, during manufacturing processes, maintenance and exploitation tasks. In engineering, the product lifecycle management (PLM) refers to the management of data and processes across the entire lifecycle of a component or a product. PLM requires a high level of data continuity across the supply chain, by developing data standards and dedicated computer platforms. Such computer platforms open the route for more machine learning in engineering, especially in MBE where engineers are facing very complex structured and unstructured data.

In this book we restrict our attention to spatially structured data. The detailed description of the data we worked with are available in Chap. 6.

### **References**


**Open Access** This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

### **Chapter 2 Learning Projection-Based Reduced-Order Models**

### **2.1 Motivation and Basic Assumptions**

In this chapter, we introduce the solution space for high-fidelity models based on partial differential equations and the finite element model. The manifold learning approach to model order reduction requires simulated data. Hence, learning projection-based reduced order models (ROM) has two steps: (i) an offline step for the computation of simulated data and for consecutive machine learning tasks, (ii) an online step where the reduced order model is used as a surrogate for the high fidelity model. The offline step generates a train set and a validation set of simulated data. The accuracy and the generalisation of the reduced order model is evaluated in the online step by using a test set of data forecast by the high-fidelity model. The test set aims also to check the computational speedups of the reduced-order model compare to the high-fidelity model.

Learning projection-based reduced order model makes sense only if there is a significant computational speedup at the price of an acceptable loss of accuracy in predictions. The longer the computational time of the high-fidelity model, the smaller the acceptable speed up, if we save hours or days of numerical simulations. Regarding the acceptable accuracy of reduced predictions, engineers working in industry and scientists working in laboratories do not have the same expectations. In our own experience, learning projection-based reduced order model has the capability to adapt to engineering tasks, high-fidelity models elaborated in laboratories, in terms of accuracy and computational time. This approach contributes to data continuity between physics-based knowledge developed in laboratories and practical applications for engineering tasks.

From the mathematical point of view, Céa's lemma gives an overview of manifold learning for model order reduction applied to elliptic equations. Let .V be a real Hilbert space, such that the weak form of the elliptic equation reads: find .~*<sup>u</sup>* <sup>∈</sup> <sup>V</sup>, such that: .*a*(~*u*,~v) <sup>=</sup> *<sup>L</sup>*(~v), <sup>∀</sup>~<sup>v</sup> <sup>∈</sup> <sup>V</sup>, (2.1)

$$a(\widetilde{u}, \widetilde{v}) = L(\widetilde{v}), \quad \forall \, \widetilde{v} \in \mathcal{V}, \tag{2.1}$$

© The Author(s) 2024

D. Ryckelynck et al., *Manifold Learning*, SpringerBriefs in Computer Science, https://doi.org/10.1007/978-3-031-52764-7\_2

where.*a*(·, ·) is a bilinear form, with coercivity constant.β > 0 and continuity continuity constant.*C<sup>a</sup>* > 0..*L*(·) is a linear form. Section 2.2 gives more details on weak forms. Let us .V*<sup>h</sup>* a finite dimensional subspace of . V. Here, .V*<sup>h</sup>* is an approximate solution space for the elliptic equation. The approximate solution of the elliptic equation is denoted by.*u* ∈ V*<sup>h</sup>* ⊂ V. The Galerkin projection of the elliptic equation onto the solution space reads: find.*u* ∈ V*<sup>h</sup>* ⊂ V such that:

$$a(u, v) = L(v), \quad \forall \ v \in \mathcal{V}\_h,\tag{2.2}$$

where.V*<sup>h</sup>* has been substituted for. V. Céa's lemma states that:

$$\text{bestimted for } \mathcal{V}. \text{ Céa's lemma states that:}$$

$$\|\widetilde{\boldsymbol{\mu}} - \boldsymbol{\mu}\| \le \frac{C^a}{\beta} \left\|\widetilde{\boldsymbol{\mu}} - \boldsymbol{\nu}\right\|, \quad \forall \, \boldsymbol{\nu} \in \mathcal{V}\_h,\tag{2.3}$$

where .|·| is a norm in . V. The related scalar product is denoted by .⟨·⟩. In Finite Element models,.V*<sup>h</sup>* is the span of finite element shape functions. But Céa's lemma holds for all finite-dimensional subspace of. V. The closer. ~*<sup>u</sup>* to the solution space.V*h*, the smaller the right-hand term of Céa's lemma (2.3), and therefore one can expect a better prediction. *<sup>u</sup>* in Eq. (2.2), although we do not know. ~*u*.

In few words, a reduced-order model is obtained by introducing a smaller solution space .V*<sup>n</sup>* ⊂ V*<sup>h</sup>* of smaller dimension .*n* < *N*. A projection-based reduced order model can be achieved by using the Galerkin projection (2.2), where.V*<sup>n</sup>* is substituted for .V*h*. The conclusion of Céa's lemma holds again. The closer . ~*<sup>u</sup>* to the solution space.V*n*, the better the prediction. *<sup>u</sup>* <sup>∈</sup> <sup>V</sup>*n*. Manifold learning comes into play when we are given a set of predictions.(*u*(*i*) )*<sup>i</sup>*=1,...,*<sup>m</sup>* in a common ambient space.V*h*, related to a given finite element mesh. The basic assumptions in manifold learning for model order reduction are:


As explained above, when the latent space is a vector subspace .V*n*, both Galerkin projection and Céa's lemma hold, but simulation speedup may not be achieved. The study of more complex situations is the purpose of this Chapter. An estimation of computational complexity of projection-based reduced order model is proposed in Sect. 2.3.6.

### **Remarks**:


### **2.2 High-Fidelity Model (HFM)**

.

~

Consider an abstract partial differential equation in a domain. Ω, with a.μ-variability: .D(~*u*; *<sup>ξ</sup>* ,*μ*) <sup>=</sup> <sup>0</sup>, *<sup>ξ</sup>* <sup>∈</sup> Ω, <sup>~</sup>*<sup>u</sup>* <sup>∈</sup> <sup>V</sup>. (2.4)

$$\mathcal{D}(\widetilde{\mu}; \xi, \mu) = 0, \quad \xi \in \Omega, \ \widetilde{\mu} \in \mathcal{V}. \tag{2.4}$$

The weak form of this partial differential equation reads: find.~*<sup>u</sup>* <sup>∈</sup> <sup>V</sup> such that

$$\text{partial differential equation reads: find } \widetilde{u} \in \mathcal{V} \text{ such that}$$

$$\int\_{\Omega} \widetilde{v} \, \mathcal{D}(\widetilde{u}; \boldsymbol{\xi}, \boldsymbol{\mu}) \, d\boldsymbol{\xi} = 0, \quad \forall \, \widetilde{v} \in \mathcal{V}. \tag{2.5}$$

As an illustration, the concepts of this chapter are illustrated on a nonlinear structural mechanics problem, for which details on the high-fidelity model are provided in this Section. For an example with another physics, we refer to the authors' work [ 21], where a nonlinear transient thermal problem is considered.

The mechanical structure occupies the domain . Ω, whose boundary .∂Ω is partitioned as.∂Ω = ∂Ω*<sup>D</sup>* ∪ ∂Ω*<sup>N</sup>* such that.∂Ω<sup>o</sup> *<sup>D</sup>* ∩ ∂Ω<sup>o</sup> *<sup>N</sup>* = ∅, see Fig. 2.1.

The structure is subjected to a quasi-static time-dependent loading, composed of an homogeneous Dirichlet boundary conditions on .∂Ω*<sup>D</sup>* and Neumann boundary conditions on .∂Ω*<sup>N</sup>* in the form of a prescribed traction .*TN* , as well as a volumic force. *f* . The setting depends on some variability. μ, which can be a parameter vector, or represent some nonparametrized variability. The evolution of the displacement .~*u<sup>μ</sup>*(*<sup>ξ</sup>* , *<sup>t</sup>*) over.(*<sup>ξ</sup>* , *<sup>t</sup>*) <sup>∈</sup> <sup>Ω</sup> × [0, *<sup>T</sup>* ] is the solution of the following equations:

$$\begin{aligned} \partial \Omega\_{D} \\\\ \epsilon(\tilde{u}^{\mu}) &= \frac{1}{2} \left( \nabla \tilde{u}^{\mu} + (\nabla \tilde{u}^{\mu})^{T} \right), & \text{in } \Omega \times [0, \mathbb{T}], \quad \text{(compatibility equation)} \\\\ \text{div} \left( \sigma\_{c}^{\mu} \right) + f^{\mu} &= 0, & \text{in } \Omega \times [0, \mathbb{T}], \quad \text{(equilibrium equation)} \\\\ \sigma\_{c}^{\mu} &= \sigma\_{c} (\epsilon(\tilde{u}^{\mu}), \mathbf{y}^{\mu}), & \text{in } \Omega \times [0, \mathbb{T}], \quad \text{(conitivity law)} \\\\ \tilde{u}^{\mu} &= 0, & \text{in } \partial \Omega\_{\mathbb{D}} \times [0, \mathbb{T}], \quad \text{(prescription law)} \\\\ \sigma\_{c}^{\mu} \cdot n\_{\partial \Omega} &= T\_{N}^{\mu}, & \text{in } \partial \Omega\_{\mathbb{D}} \times [0, \mathbb{T}], \quad \text{(prescription law)} \\\\ \tilde{u}^{\mu} &= 0, \, \text{y}^{\mu} = 0, & \text{in } \Omega \text{ at } t = 0, \quad \text{(initial condition)} \end{aligned}$$

)

(2.6)

where . *∈* is the linear strain tensor, .*σ <sup>μ</sup> <sup>c</sup>* is the Cauchy stress tensor, .*y<sup>μ</sup>* denotes the internal variables of the constitutive law and .*n*∂Ω is the outward normal vector on .∂Ω. We precise that the evolution of the internal variables .*y<sup>μ</sup>* is updated when the constitutive law is solved.

Define .*H*<sup>1</sup> <sup>0</sup> (Ω) = {<sup>v</sup> <sup>∈</sup> *<sup>L</sup>*<sup>2</sup>(Ω)<sup>|</sup> ∂v ∂ξ*<sup>i</sup>* <sup>∈</sup> *<sup>L</sup>*<sup>2</sup>(Ω), <sup>1</sup> <sup>≤</sup> *<sup>i</sup>* <sup>≤</sup> 3 and <sup>v</sup>|∂Ω*<sup>D</sup>* <sup>=</sup> <sup>0</sup>}. Denote .{ϕ*i*}1≤*i*≤*<sup>N</sup>* <sup>∈</sup> <sup>R</sup>*N*×*<sup>N</sup>* , a finite element basis whose span, denoted .V*h*, constitutes an approximation of.*H*<sup>1</sup> <sup>0</sup> (Ω)3;.*N* is the number of finite element basis functions, hence the number of degrees of freedom of the discretized prediction. A discretized weak formulation reads: find.*u<sup>μ</sup>* ∈ V*<sup>h</sup>* such that for all.v ∈ V*h*, {{{

$$\int\_{\Omega} \sigma\_{\varepsilon}(\epsilon(u^{\mu}), \mathbf{y}) : \epsilon(v) = \int\_{\Omega} f^{\mu} \cdot v + \int\_{\partial \Omega\_{N}} T\_{N}^{\mu} \cdot v,\tag{2.7}$$

that we denote for concision.**F***<sup>μ</sup>* (**u***<sup>μ</sup>*) <sup>=</sup> 0, where.**u***<sup>μ</sup>* is the vector of.*<sup>N</sup>* coordinates for .*u<sup>μ</sup>* ∈ V*h*. A Newton algorithm can be used to solve this nonlinear global equilibrium problem at each time step: (**<sup>u</sup>**μ,*<sup>k</sup>* ) (**u***<sup>μ</sup>*,*k*+<sup>1</sup> <sup>−</sup> **<sup>u</sup>***<sup>μ</sup>*,*<sup>k</sup>* ) = −**F***μ* ( **u***μ*,*k* )

$$
\frac{D\mathbf{F}^{\mu}}{Du} \left(\mathbf{u}^{\mu,k}\right) \left(\mathbf{u}^{\mu,k+1} - \mathbf{u}^{\mu,k}\right) = -\mathbf{F}^{\mu} \left(\mathbf{u}^{\mu,k}\right), \tag{2.8}
$$

$$
\frac{\mu}{\mu} \left(\mathbf{u}^{k}\right)\_{\shortdot} = \int\_{\alpha} \epsilon \left(\boldsymbol{\varphi}\_{j}\right) : \mathcal{K} \left(\boldsymbol{\epsilon}(\boldsymbol{u}^{\mu,k}), \mathbf{y}^{\mu}\right) : \boldsymbol{\epsilon} \left(\boldsymbol{\varphi}\_{i}\right), \tag{2.9}
$$

where

$$\frac{D\mathbf{F}^{\mu}}{Du}\left(\mathbf{u}^{k}\right)\_{ij} = \int\_{\Omega} \boldsymbol{\epsilon}\left(\boldsymbol{\varphi}\_{j}\right) : \mathcal{K}\left(\boldsymbol{\epsilon}(\boldsymbol{u}^{\mu,k}), \mathbf{y}^{\mu}\right) : \boldsymbol{\epsilon}\left(\boldsymbol{\varphi}\_{i}\right),\tag{2.9}$$

and

(

.

.

.

)

.

**Fig. 2.2** Linear manifold learning

\*\*p.2.2\*\*. Linear manifold learning

$$\mathbf{F}^{\mu} \left( \mathbf{u}^{\mu,k} \right)\_{i} = \int\_{\Omega} \sigma\_{c} \left( \epsilon(\boldsymbol{u}^{\mu,k}), \mathbf{y}^{\mu} \right) : \epsilon \left( \boldsymbol{\varrho}\_{i} \right) - \int\_{\Omega} f^{\mu} \cdot \boldsymbol{\varrho}\_{i} - \int\_{\partial \Omega\_{N}} T\_{N}^{\mu} \cdot \boldsymbol{\varrho}\_{i}. \tag{2.10}$$

In the two relations above,  $\mathcal{K} \left( \epsilon(\boldsymbol{u}^{\mu,k}), \mathbf{y}^{\mu} \right)$  is the local tangent operator,  $\boldsymbol{\mu}^{\mu,k} \in \mathcal{V}$ .

is the local tangent operator,. *uμ*,*<sup>k</sup>* ∈ V is the k-th iteration of the discretized displacement field at the current time-step, and <sup>∈</sup> <sup>R</sup>*<sup>N</sup>* is such that.*uμ*,*<sup>k</sup>* <sup>=</sup> <sup>∑</sup> ()(

.**u***μ*,*<sup>k</sup>* = *uμ*,*<sup>k</sup> i* 1≤*i*≤*N N i*=1 *uμ*,*<sup>k</sup> <sup>i</sup>* <sup>ϕ</sup>*<sup>i</sup>* . In particular,. *<sup>f</sup> <sup>μ</sup>*,.*<sup>T</sup>* <sup>μ</sup> *<sup>N</sup>* ,. *uμ*,*<sup>k</sup>* ), *yμ*) (

and.*y<sup>μ</sup>* are known and enforce the time-dependence of the solution. Depending on the constitutive law, the computation of the functions. *uμ*,*<sup>k</sup>* , *y* |→ *σ<sup>c</sup> ∈*(*uμ*,*<sup>k</sup>* and . *uμ*,*<sup>k</sup>* , *yμ*) |→ K ( *∈*(*uμ*,*<sup>k</sup>* ), *yμ*) can require the resolution of a complex differentialalgebraic system of equations.

### **2.3 Linear Manifold Learning for Projection-Based Reduced-Order Modeling**

We start by explaining the online phase. Since we want to construct and solve the reduced-order model (ROM) in the most efficient way, the offline phase is dedicated to precompute as many steps as possible, under the considered variability.

Linear manifold learning means that the solution manifold is approximated by a vector subspace of the ambient solution space, as illustrated in Fig. 2.2.

### *2.3.1 Approaches Preceding the Use of Machine Learning*

In structural mechanics, normal modes have been introduced for the analysis of vibrations in structures. When considering free vibrations, without external force, the solution of linear hyperbolic equation is sought by using the separation of space and time variables: .**u***<sup>N</sup>* (**x**, *<sup>t</sup>*) <sup>=</sup> *<sup>ψ</sup>*(**x**) *<sup>u</sup>*(*t*), where .**<sup>x</sup>** <sup>∈</sup> <sup>Ω</sup> is the space variable and . *<sup>t</sup>* <sup>∈</sup> <sup>R</sup> the time variable. The hyperbolic equation of free vibration reads: find. *ψ*(**x**) ∈ V*<sup>N</sup>* and. *<sup>u</sup>*(*t*) <sup>∈</sup> <sup>R</sup> such that.**u***<sup>N</sup>* (**x**, *<sup>t</sup>*) <sup>=</sup> *<sup>ψ</sup>*(**x**) *<sup>u</sup>*(*t*) and

$$
\langle \rho \: \ddot{\mathbf{u}}\_N, \mathbf{v}\_N \rangle + a(\mathbf{u}\_N, \mathbf{v}\_N) = 0, \quad \forall \, \mathbf{v}\_N \in \mathcal{V}\_N. \tag{2.11}
$$

It follows that . *<sup>u</sup>*(*t*) is an harmonic function of frequency denoted by . *<sup>f</sup>* and .*<sup>ψ</sup>* is the eigenvector related to the eigenvalue.λ = (2 π *f* )<sup>2</sup> of the following generalized eigenproblem: find.*ψ* ∈ V*<sup>N</sup>* and. λ such that

$$a(\boldsymbol{\Psi}, \mathbf{v}\_N) - \lambda \,\,\langle \boldsymbol{\rho} \,\, \boldsymbol{\Psi}, \mathbf{v}\_N \rangle = 0, \quad \forall \, \mathbf{v}\_N \in \mathcal{V}\_N,\tag{2.12}$$

where the rank of this system of equation is supposed to be.*N* − 1 in order to find non zero eigenmodes. This eigenproblem admits.*N* orthogonal normal modes that span .V*<sup>N</sup>* . Therefore a selection of. *n* normal modes span a reduced subspace of dimension . *n*. In the beginning of the . 21st century, model reduction using variable separation scheme in partial differential equations has been extended to an arbitrary number of variable by using low-rank approximations such as the Proper Generalized Decomposition [ 5]. Eigenmodes are global functions in contrast to finite element shape functions that have a local support. For dynamical problems involving nonlinear contributions to the PDE, adaptive computations of reduced subspaces have been proposed by Almroth et al. [ 4] and Noor et al. [ 79], by using Rayleigh-Ritz global functions. The set of these global functions is a reduced basis of the finite-element approximation-space.

The idea of using statistics to generate a solution space for differential equations has been proposed in the seminal work of Lorenz in [ 69] (in page 31), by using empirical orthogonal functions. The Galerkin projection of PDEs on empirical modes have been first developed in [ 13], where a reduced basis is computed via the proper orthogonal decomposition [ 70] of observational data. This was the first step towards manifold learning for projection-based model order reduction, which we now present. For other presentations of this technologies, the reader can refer to [ 81, 87].

### *2.3.2 Online Phase: Galerkin Projection*

.

The reduced-order model is constructed in the form of a Galerkin method written on a Reduced-Order Basis (ROB). In the present case, it consists in assembling the physical problem in the same fashion as the HFM in Sect. 2.2, with the difference that the finite element basis.(ϕ*i*)<sup>1</sup>≤*i*≤*<sup>N</sup>* <sup>∈</sup> <sup>R</sup>*<sup>N</sup>*×*<sup>N</sup>* , is replaced by a ROB.(ψ*i*)<sup>1</sup>≤*i*≤*<sup>n</sup>* <sup>∈</sup> <sup>R</sup>*<sup>n</sup>*×*<sup>N</sup>* , with.*n* < *N*. Hence, the reduced Newton algorithm is constructed as () (*<sup>γ</sup> <sup>μ</sup>*,*k*+<sup>1</sup> <sup>−</sup> *<sup>γ</sup> <sup>μ</sup>*,*<sup>k</sup>* ) ()

$$\frac{D\mathcal{F}\_{\mu}}{D\hat{u}}\left(\hat{u}\_{\mu}^{k}\right)\left(\mathcal{Y}^{\mu,k+1} - \mathcal{Y}^{\mu,k}\right) = -\mathcal{F}\_{\mu}\left(\hat{u}\_{\mu}^{k}\right),\tag{2.13}$$

### 2.3 Linear Manifold Learning for Projection-Based … 15 (){()

where

(

)

$$
\begin{aligned}
\text{ou Lemmaing to 1 tojecauon-basau }\dots \\
\\
\frac{D\overline{\mathcal{F}}\_{\mu}}{D\hat{u}}\left(\hat{u}\_{\mu}^{k}\right)\_{ij} = \int\_{\Omega} \epsilon\left(\psi\_{j}\right) : \mathcal{K}\left(\epsilon(\hat{u}\_{\mu}^{k}), \mathbf{y}\_{\mu}\right) : \epsilon\left(\psi\_{i}\right)\end{aligned}
\tag{2.14}
$$

and

$$
\hat{D}\hat{\boldsymbol{\mu}}^{k} \stackrel{\leftarrow}{\quad} \hat{\boldsymbol{\mu}}^{k} \boldsymbol{\mu}\_{ij} \quad \boldsymbol{J}\_{\Omega} \stackrel{\leftarrow}{\quad} \hat{\boldsymbol{\epsilon}}^{k,N} \quad \boldsymbol{\epsilon}^{k} \boldsymbol{\epsilon}^{k} \boldsymbol{\mu}^{k} \boldsymbol{\epsilon}^{k} \boldsymbol{\mu}^{\prime} \quad \boldsymbol{\epsilon}^{\prime} \boldsymbol{\epsilon}^{k}
$$

$$
\mathcal{F}\_{\mu} \left(\hat{\boldsymbol{\mu}}^{k}\_{\mu}\right)\_{i} = \int\_{\Omega} \sigma\left(\boldsymbol{\epsilon}(\hat{\boldsymbol{\mu}}^{k}\_{\mu}), \mathbf{y}\right) : \boldsymbol{\epsilon}\left(\boldsymbol{\psi}\_{i}\right) - \int\_{\Omega} f\_{\mu} \cdot \boldsymbol{\psi}\_{i} - \int\_{\partial\Omega\_{N}} T\_{\mu,N} \cdot \boldsymbol{\psi}\_{i}. \tag{2.15}
$$

In the two relations above, .*u*ˆ *<sup>k</sup>* <sup>μ</sup> ∈ Vˆ := Span (ψ*i*)<sup>1</sup>≤*i*≤*<sup>n</sup>* is the .*k*-th iteration of the reduced displacement field for the current time-step and . *γ <sup>μ</sup>*,*<sup>k</sup>* = γ *<sup>μ</sup>*,*<sup>k</sup> i* <sup>1</sup>≤*i*≤*<sup>n</sup>* <sup>∈</sup> <sup>R</sup>*<sup>n</sup>* is such that μ = ∑*n*

$$
\hat{\mu}^k\_{\mu} = \sum\_{i=1}^n \nu\_i^{\mu,k} \psi\_i. \tag{2.16}
$$

)

Notice that the use of the Galerkin method is made possible by the linearity of the tangent problem (2.13) and the choice of a linear dimensionality reduction technique in (2.16).

The online stage is called *efficient* if the reduced problem can be constructed and solved in computational complexity independent of . *N*. When the variability . μ is parametrized, efficiency is possible by precomputing various terms. With nonparametrized variability, depending on its nature, some assembling task with a linear complexity in .*N* may be required at the beginning of the online stage (for instance for a boundary condition). All these scenarios are handled by genericROM, the ROM library developped at Safran and presented in Sect. 4.2.

The offline phase contains three steps, for which we present below the methodological choices made in genericROM.

### *2.3.3 Offline Phase*

### **2.3.3.1 Data Generation**

This step corresponds to the generation of the snapshots by solving the high-fidelity model. In parametric contexts, the simplest workflows consist in choosing parameter values a priori, following Design of Experiments (DoE, see [ 50] for a recent technique), and computing the corresponding snapshots by solving the high-fidelity model.

### **2.3.3.2 Data Compression: Dimensionality Reduction**

This step corresponds to the generation of the ROB.(ψ*i*)1≤*i*≤*<sup>n</sup>*. One of the most classical and simple method is the snapshot Proper Orthogonal Decomposition (POD) [ 25, 86], detailed below: 2. Compute the correlation matrix.*Ci*,*<sup>j</sup>* <sup>=</sup> {


$$\text{4. Compute the reduced order basis } \psi\_i(\mathbf{x}) = \frac{1}{\sqrt{\lambda\_i}} \sum\_{j=1}^{N\_s} u\_j(\mathbf{x}) \xi\_{i,j}, 1 \le i \le n.$$

The advantages of the snapshot-POD are a reasonable computational complexity when the number of degree of freedom of the high-fidelity model are much larger than the number of snapshots, and the fact that this algorithm can be easily parallelized.

Variants can be used, for instance in the Reduced Basis [ 84] and the PODgreedy [ 41] methods, with respectively an orthonormalization of the computed highfidelity snapshots and the incremental Singular Value Decomposition (SVD) [ 15].

### **2.3.3.3 Operator Compression**

A ROM is called online-efficient if in the online stage, the reduced problems can be constructed and solved in computational complexity independent of. *N*. The operator compression step consists in additional treatments required for the efficiency of the online stage, by pre-processing some computationally demanding integration tasks over the high-fidelity domain .Ω and .∂Ω*<sup>N</sup>* . We notice that without any additional treatment, the numerical integration involved in the assembling of Eq. (2.13) strongly limits in practice the efficiency of the ROM: no speedup with respect to the highfidelity model can be obtained in practice. The complexity of such an additional treatment depends on the type of parameter-dependence of the problem. This step is actually needed for all classes of problems reduced by projection-based methods.

Consider the simplest case: a linear problem with an affine dependence in the parameter . μ, for instance .*A*μ**u** = **c**, where .*A*<sup>μ</sup> = *A*<sup>0</sup> + μ*A*1. Denote **V**, the matrix whose columns are the vectors of the ROB evaluated at the high-fidelity degrees of freedom. The obtained ROM writes .**V***<sup>T</sup> A*μ**Vu**ˆ = **V***<sup>T</sup>* **c**: it is not assembled in the online phase, but rather the matrices .**V***<sup>T</sup> A*0**V** and .**V***<sup>T</sup> A*1**V** and the vector .**V***<sup>T</sup>* **c** are precomputed in the offline stage so that the reduced problem is constructed without approximation and efficiently by summing two small matrices. The operator compression step consists in this case in the construction of.**V***<sup>T</sup> A*0**V** and.**V***<sup>T</sup> A*1**V** in the offline stage.

Actually, there exist linear problems for which the operation compression step require an additional approximation, and nonlinear problems that can be carried-

out exactly. In the first case, consider .*A*μ**u** = **c** with . *Ai j* = Ω ∇ ( *g*(*x*, μ)ϕ*j*(*x*) · <sup>∇</sup>ϕ*i*(*x*) and .*ci* <sup>=</sup> { <sup>Ω</sup> *f* (*x*)ϕ*i*(*x*), where . *u* is the unknown, . *f* a known loading and .*g*(*x*, μ) a known function with variables . *x* and .μ that cannot be separated: the previous precomputation of reduced matrices cannot be applied, and a treatment is required to, for example, approximately separate the dependencies in . *x* and . μ of . *g* as .*g*(*x*, μ) <sup>≈</sup> <sup>∑</sup>*<sup>d</sup> <sup>k</sup>*=<sup>1</sup> *<sup>g</sup><sup>a</sup> <sup>k</sup>* (*x*)*g<sup>b</sup> <sup>k</sup>* (μ). Then, .**V***<sup>T</sup> <sup>A</sup>*μ**<sup>V</sup>** <sup>≈</sup> <sup>∑</sup>*<sup>d</sup> <sup>k</sup>*=<sup>1</sup> *<sup>g</sup><sup>b</sup> <sup>k</sup>* (μ)*Ak* where . (*Ak* )*i j* <sup>=</sup> { Ω ∇ ( *ga <sup>k</sup>* (*x*)ϕ*j*(*x*) ) · ∇ϕ*i*(*x*), so that the efficiency of the online stage is recovered; the Empirical Interpolation Method has been proposed in [14, 71] for this purpose. A case of linear problem in harmonic aeroacoustics with nonaffine dependence with respect to the frequency is available in [ 20]. Conversely, nonlinearities can be handled without approximation in some cases, for instance, the advection term in fluid dynamics can be precomputed in the form of an order-3 tensor:. { Ω ψ*<sup>i</sup>* · ( ψ*j* · ∇) ψ*<sup>k</sup>* ,.1 ≤ *i*, *j*, *k* ≤ *n*; see [ 3] for the reduction of the nonlinear Navier-Stokes equations with an exact operator compression step. Other examples are found in structural dynamics with geometric nonlinearities, where order-2 and -3 tensors can also be precomputed, see [ 58, Sect. 3.2] and [ 73].

When additional approximations are required, the methods proposed for the operator compression step are call "hyper-reduction" in the literature. This term was coined by the seminal method proposed in [ 88] in 2005, but has been extended to refer to all the methods proposing a such second reduction stage. Hyper-reduction methods include the Empirical Interpolation Method (EIM, [ 14]), the Missing Point Estimation (MPE, [ 12]), the Best Point Interpolation Method (BPIM, [ 78]), the Discrete Empirical Interpolation Method (DEIM, [ 23]), the Gauss-Newton with Approximated Tensors (GNAT, [ 17]), the Energy-Conserving Sampling and Weighting (ECSW, [ 36]), the Empirical Cubature Method (ECM, [ 45]), and the Linear Program Empirical Quadrature Procedure (LPEQP, [100]). The reader can find an algorithmic comparison of the Hyper-Reduction and the Discrete Empirical Interpolation Method for a nonlinear thermal problem in [ 39]. A particular focus is given on hyper-reduction techniques via oblique projection and empirical cubature in the following sections.

### *2.3.4 Hyper-Reduction via a Reduced Integration Domain*

Hyper-reduction via a reduced integration domain has been proposed in [ 88]. It requires a train set of displacement predictions, so that a reduced approximation vector space can be trained. The finite elements simulations that generate the train set of displacement fields are also predicting stresses fields .*σ* and internal variables. **y**. Hence, additional reduced bases can be trained for these variables, by using these simulation results [ 89]. Heuristically, we found it more accurate to include a reduced basis for the stresses. *σ*, as an additional reduced basis for this hyper-reduction

{

scheme. Such a reduced basis is also very convenient for error estimation [ 91]. The finite element shape functions for displacement fields are denoted by .(*ϕi*)*<sup>i</sup>*=1,...,*<sup>N</sup>* . For stress fields, we also need to introduce a related finite element representation. We denote by .(*ϕ*<sup>σ</sup> *<sup>i</sup>* )*<sup>i</sup>*=1,...,*N*<sup>σ</sup> the dedicated shape functions. In the linear framework of manifold learning, we assume that the same finite element mesh is used for the target simulation, in the online step or for the test set of data, and all simulations used to generate the train set of data.

In practice, the implementation of the hyper-reduction follows the manifold learning step that trains reduced bases for displacements and stresses. We recall that they are respectively denoted by .**<sup>V</sup>** <sup>∈</sup> <sup>R</sup>*N*×*<sup>n</sup>* and .**V**<sup>σ</sup> <sup>∈</sup> <sup>R</sup>*N*<sup>σ</sup> <sup>×</sup>*n*<sup>σ</sup> in their matrix form, and .(*ψ<sup>k</sup>* )*<sup>k</sup>*=1,...,*<sup>n</sup>* .(*ψ*<sup>σ</sup> *<sup>k</sup>* )*<sup>k</sup>*=1,...,*n*<sup>σ</sup> in their continuous form: .*ψk* (**x**) = ∑

$$\Psi\_k(\mathbf{x}) = \sum\_{i=1}^N \varphi\_i(\mathbf{x}) V\_{ik}, \quad \forall \mathbf{x} \in \Omega, \ k = 1, \ldots, n,\tag{2.17}$$

$$\boldsymbol{\Psi}\_{k}^{\sigma}(\mathbf{x}) = \sum\_{i=1}^{N^{\sigma}} \boldsymbol{\varphi}\_{i}^{\sigma}(\mathbf{x}) V\_{ik}^{\sigma}, \quad \forall \mathbf{x} \in \mathfrak{Q}, \ k = 1, \ldots, n^{\sigma}. \tag{2.18}$$

The reduced displacement reads:

$$
\widehat{\mathbf{u}} = \mathbf{u}
$$

$$
\text{hent reads:}
$$

$$
\widehat{\mathbf{u}}(\mathbf{x}) = \sum\_{k=1}^{n} \boldsymbol{\Psi}\_{k}(\mathbf{x}) \ \boldsymbol{\eta}\_{k}, \ \forall \mathbf{x} \in \Omega.
\tag{2.19}
$$

The hyper-reduction method proposed in [88] aims at computing reduced coordinates .(γ*<sup>k</sup>* )*<sup>k</sup>*=1,...,*<sup>n</sup>* introduced in Eq. (2.19), by projecting the equilibrium equation on. **V**, via a restriction of the domain .Ω to a Reduced Integration Domain (RID) denoted by .Ω*R*. By following the empirical interpolation method [ 14], interpolation points are computed for column vectors in .**V** and .**V**<sup>σ</sup> separately [ 48]. The set of respective interpolation points are denoted by.P*<sup>u</sup>* and.P<sup>σ</sup> . We follow Algorithm 2.1, proposed for the Discrete Empirical Interpolation Method [ 23].

The RID.Ω*<sup>R</sup>* is such that it contains the interpolation points related to. **V** and.**V**<sup>σ</sup> . For engineering applications, the RID can also include a zone of interest in .Ω that is user-defined, by using a subset of finite elements. This zone of interest is denoted by.Ω*Z I* ⊂ Ω. By construction, for contact-free problems, the RID is the following:

$$
\Omega\_R = \Omega\_{ZI} \cup\_{i \in \mathcal{P}^w} \text{supp}(\mathfrak{\mathfrak{p}}\_i) \cup\_{i \in \mathcal{P}^w} \text{supp}(\mathfrak{\mathfrak{q}}\_i^{\sigma}), \tag{2.20}
$$

where.supp( *f* ) ∈ Ω is the support of function. *f* . In practice,.Ω*<sup>R</sup>* has its own finiteelement mesh. It is a reduced mesh involving much less elements than the original finite element mesh of. Ω. It can be enlarge by adding a layer of connected element in the original mesh. Integration of constitutive equations are performed on this reduced mesh, without any intrusive operation on the original finite element solver. A similar hyper-reduction scheme has been developed in [ 38] for contact problems.

{

.

{

**Algorithm 2.1:** Interpolation points of the Discrete Empirical Interpolation Method (DEIM) [ 23]

**Input :** reduced basis vectors **V**[:, *<sup>k</sup>*] ∈ <sup>R</sup>*N* , *k* <sup>=</sup> <sup>1</sup>,... *<sup>M</sup>* **Output:** interpolation point index set *I* := {*i*1,...,*iM* } set *I*0 := ∅ ; // initialization **2 for** *l* = 1 ..., *M* **do U***l*−1 ← **V**[:, 1 : (*l* − 1)] ; // truncated basis **<sup>A</sup>**<sup>←</sup> (**U***l*−1[*Il*−1, :]*T* **U***l*−1[*Il*−1, :])−1**U***l*−1[*Il*−1, :]*T* ; // projector **q***l* ← **V**[:,*l*] − **U***l*−1**AV**[*Il*−1,*l*] ; // interpolation residual *il* ← arg max*i*∈{1,...,*N* } |(**q***l* )*i*| ; // maximum of residual *Il* := *Il*−1 ∪ {*il*} ; // extend interpolation points **8 end**  set *I* := *IM* .

Once the RID is obtained, a set of test functions is set up in order to restrain the balance equations to.Ω*R*. They are denoted by.*ψ R j* : {

$$\mathcal{P} = \left\{ i \in \{ 1, \ldots, N \}, \int\_{\Omega \backslash \Omega\_{\mathbb{R}}} \left( \mathfrak{p}\_{i} \right)^{2} d\Omega = 0 \right\},\tag{2.21}$$

}

$$\mathcal{P} = \left\{ i \in \{ 1, \ldots, N \} , \int\_{\Omega \backslash \Omega\_{\mathbb{R}}} (\mathfrak{p}\_{i})^{2} \, d\Omega = 0 \right\},\tag{2.21}$$

$$\Psi\_{R, j}(\mathbf{x}) = \sum\_{i \in \mathcal{P}} \mathfrak{p}\_{i}(\mathbf{x}) \, V\_{ij}, \quad \forall \mathbf{x} \in \Omega, \,\, j = 1, \ldots, n,\tag{2.22}$$

where .P is the set of all degrees of freedom in .Ω*<sup>R</sup>* excepted those belonging to the interface between .Ω*<sup>R</sup>* and its counterpart. This interface is denoted by .I*R*. As explained in [ 94], the test functions are null on the interface.I*R*, as if Dirichlet boundary conditions were imposed to the RID. On this interface, displacements follow the shape of the modes .*ψ<sup>k</sup>* according to Eq. (2.19). The hyper-reduction method gives access to reduced coordinates.(γ*<sup>k</sup>* )*<sup>k</sup>*=1,...,*<sup>n</sup>* that fulfill the following balance equations, for contactless problems: . **u**(**x**) = ∑*n*

$$
\widehat{\mathbf{u}}(\mathbf{x}) = \sum\_{k=1}^{n} \boldsymbol{\Psi}\_{k}(\mathbf{x}) \, \boldsymbol{\uprho}\_{k} \, , \, \forall \mathbf{x} \in \Omega\_{R} \tag{2.23}
$$

$$
\int\_{\Omega\_{R}} \mathbf{e}(\boldsymbol{\upPsi}\_{R,j}) \, : \, \boldsymbol{\sigma}(\mathbf{e}(\widehat{\mathbf{u}})) \, d\Omega \tag{2.24}
$$

$$\int\_{\Omega\_{\mathcal{R}}} \mathfrak{s}(\boldsymbol{\Psi}\_{\mathcal{R}\_{\mathcal{J}}}) \; : \; \mathfrak{o}(\mathfrak{e}(\widehat{\mathfrak{u}})) \; d\Omega \tag{2.24}$$

$$-\int\_{\Omega\_R} \Psi\_{R,j} \, f\_{\mu} \, d\Omega - \int\_{\partial\Omega\_R \cap \partial\Omega\_N} \Psi\_{R,j} \, T\_{\mu} \, N \, dS = 0. \tag{2.25}$$

$$\forall j = 1, \ldots, n \tag{2.26}$$

The matrix form of the hyper-reduced balance equations reads: find .*<sup>γ</sup>* <sup>∈</sup> <sup>R</sup>*<sup>n</sup>* such that

$$
\widehat{\mathbf{u}}(\mathbf{x}) = \sum\_{i=1}^{N} \boldsymbol{\varphi}\_{i}(\mathbf{x}) \,\widehat{q}\_{i}(\mathbf{x}) \,\widehat{q}\_{i}(\mathbf{x}) \,\widehat{q}\_{i}(\mathbf{x}) \,\widehat{q}\_{i}(\mathbf{x}) \,\widehat{q}\_{i}(\mathbf{x}) \,\widehat{q}\_{i}(\mathbf{x}) \,\widehat{q}\_{i}(\mathbf{x}) \,\widehat{q}\_{i}(\mathbf{x}) \,\widehat{q}\_{i}(\mathbf{x}) \,\widehat{q}\_{i}(\mathbf{x}) \,\widehat{q}\_{i}(\mathbf{x}) \,\widehat{q}\_{i}(\mathbf{x}) \,\widehat{q}\_{i}(\mathbf{x}) \,\widehat{q}\_{i}(\mathbf{x}) \,\widehat{q}\_{i}(\mathbf{x}) \,\widehat{q}\_{i}(\mathbf{x}) \,\widehat{q}\_{i}(\mathbf{x}) \,\widehat{q}\_{i}(\mathbf{x}) \,\widehat{q}\_{i}(\mathbf{x}) \,\widehat{q}\_{i}(\mathbf{x}) \,\widehat{q}\_{i}(\mathbf{x}) \,\widehat{q}\_{i}(\mathbf{x}) \,\widehat{q}\_{i}(\mathbf{x}) \,\widehat{q}\_{i}(\mathbf{x}) \,\widehat{q}\_{i}(\mathbf{x}) \,\widehat{q}\_{i}(\mathbf{x}) \,\widehat{q}\_{i}(\mathbf{x}) \,\widehat{q}\_{i}(\mathbf{x}) \,\widehat{q}\_{i}(\mathbf{x}) \,\widehat{q}\_{i}(\mathbf{x}) \,\widehat{q}\_{i}(\mathbf{x}) \,\widehat{q}\_{i}(\mathbf{x}) \,\widehat{q}\_{i}(\mathbf{x}) \,\widehat{q}\_{i}(\mathbf{x}) \,\widehat{q}\_{i}(\mathbf{x}) \,\widehat{q}\_{i}(\mathbf{x}) \,\widehat{q}\_{i}(\mathbf{x}) \,\$$

$$
\widehat{\mathbf{q}} = \mathbf{V}\,\mathbf{y},\tag{2.28}
$$

$$\mathcal{F}^{HR}(\mathbf{y}) \coloneqq \mathbf{V}[\mathcal{P}, \mathbf{:}]^T \mathcal{F}(\mathbf{V} \mathbf{y})[\mathcal{P}],\tag{2.29}$$

$$\mathcal{F}^{HR}(\mathfrak{y}) = 0,\tag{2.30}$$

where .**V**[P, :] denotes a row restriction of matrix .**V** to indices in . P. The reduced Newton-Raphson step reads:

$$\begin{aligned} \text{row restriction of matrix } \mathbf{V} \text{ to indices in } \mathcal{P}. & \text{The reduced} \\\\ \widehat{\mathbf{u}}^{k}(\mathbf{x}) &= \sum\_{i=1}^{N} \boldsymbol{\varphi}\_{i}(\mathbf{x}) \, \widehat{q}\_{i}^{k}, \,\forall \mathbf{x} \in \Omega\_{R}, \\\ \widehat{\mathbf{q}}^{k} &= \mathbf{V} \, \mathbf{y}^{k}, \end{aligned} \tag{2.31}$$

$$\begin{aligned} \mathbf{\hat{v}}^{k} &:= \sum\_{i=1}^{I} \mathbf{\hat{v}}^{k} [\mathbf{\hat{v}}^{k}, \mathbf{\hat{u}}], \quad \mathbf{\hat{u}} &:= \mathbf{\hat{u}}^{k}, \mathbf{\hat{v}}^{k} \\ \mathbf{\hat{q}}^{k} &:= \mathbf{V} \mathbf{\hat{y}}^{k}, \quad \mathbf{\hat{z}}^{k} := \mathbf{V} [\mathcal{P}, \mathbf{:}]^{T} \frac{D \mathcal{F}\_{\mu}}{D \mathbf{\hat{x}}} (\mathbf{\hat{u}}^{k-1}) [\mathcal{P}, \mathbf{:}] \mathbf{V}, \quad \mathbf{\hat{z}}^{k} \end{aligned} \tag{2.32}$$

$$\begin{aligned} \widehat{\mathbf{q}}^k &= \mathbf{V} \mathbf{y}^k, \\ \mathbf{K}^{HR} := \mathbf{V}[\mathcal{P}, :]^T \frac{D\mathcal{F}\_\mu}{D\widehat{\mathbf{u}}}(\widehat{\mathbf{u}}^{k-1})[\mathcal{P}, :] \, \mathbf{V}, \\ \mathbf{y}^{k-1}) &= -\mathbf{V}[\mathcal{P}, :]^T \mathcal{F}(\widehat{\mathbf{u}}^{k-1})[\mathcal{P}], \end{aligned} \tag{2.33}$$

$$\mathbf{K}^{HR} \left( \mathbf{y}^{k} - \mathbf{y}^{k-1} \right) = -\mathbf{V}[\mathcal{P}, \mathbf{:}]^T \mathcal{F}(\widehat{\mathbf{u}}^{k-1})[\mathcal{P}],\tag{2.34}$$

where the reduced stiffness matrix.**K***H R* is computed by using solely the elements of the RID.Ω*R*. We assume that the matrix.**K***H R* is full rank. This assumption is always checked in numerical solutions of hyper-reduced equations. Rank deficiency may appear when the RID construction do not account for the contribution of a reduced basis dedicated to stresses.

Once the RID is represented as a finite element mesh, this hyper-reduction scheme is intrusive solely for the linear solver involved in the Newton-Raphson step and its related convergence criterion. Nevertheless, the mesh of the RID has to include labels for the set. P or its counterpart.I*R*. This counterpart is the set of degrees of freedom connected to elements of the original mesh that are not in the reduced mesh.

### **Remarks**:


.

∑

• Reduced order models not only save computational time, they save computational resources including energy consumption savings as explained in [ 90] and memory footprint [ 46].

**Property 1**: In linear elasticity, if .**K***H R* is full rank, the hyper-reduced balance equations are equivalent to an oblique projection of the finite element prediction .**<sup>q</sup>** <sup>∈</sup> <sup>R</sup>*<sup>N</sup>* :

$$\mathbf{II}^T \coloneqq \mathbf{V}[\mathcal{P}, \mathbf{:}]^T \mathbf{K}[\mathcal{P}, \mathbf{:}],\tag{2.35}$$

$$\begin{aligned} \boldsymbol{\Pi}^{T} &:= \mathbf{V}[\mathcal{P}, \mathbf{:}]^{T} \, \mathbf{K}[\mathcal{P}, \mathbf{:}], \\ \widehat{\mathbf{q}} &= \mathbf{V} \, (\boldsymbol{\Pi}^{T} \, \mathbf{V})^{-1} \boldsymbol{\Pi}^{T} \, \mathbf{q}, \\ \text{and} \quad \boldsymbol{\Pi}^{T} \, \widehat{\mathbf{q}} &= \boldsymbol{\Pi}^{T} \, \mathbf{q}, \end{aligned} \tag{2.35}$$

)

$$\mathbf{I} \text{ and } \begin{aligned} \Pi^I \ \widehat{\mathbf{q}} = \Pi^I \ \mathbf{q}, \end{aligned} \tag{2.37}$$

with .**K q** = **F**. Hence the hyper-reduced prediction of the reduced vector .*γ* is a minimizer for. *f* (*β*): ) = |*∏T* (

$$\mathbf{y}^\star \in \mathbb{R}^n, \ f(\mathbf{y}^\star) = \|\boldsymbol{\Pi}^T \left(\mathbf{V} \,\mathbf{y}^\star - \mathbf{q}\right)\|\_2^2. \tag{2.38}$$

Here.*∏* is a projector for elastic stresses in.Ω*<sup>R</sup>* according to the reduced test functions:

$$\begin{aligned} \text{projection for elastic stresses in } \Omega\_R \text{ according to the reduced test functions:}\\ \sum\_{i=1}^N \Pi\_{ik} \left( \mathbf{V} \, \mathbf{y} - \mathbf{q} \right)\_i = \int\_{\Omega\_R} \mathbf{e} \left( \boldsymbol{\Psi}\_{R,k} \right) \, : \, \left( \boldsymbol{\sigma}(\widehat{\mathbf{q}}) - \boldsymbol{\sigma}(\mathbf{q}) \right) \, d\Omega. \end{aligned} \tag{2.39}$$

The proof is straightforward. Here, .**K***H R* <sup>=</sup> *<sup>∏</sup><sup>T</sup>* **<sup>V</sup>**. The Jacobian matrix for . *<sup>f</sup>* reads .**<sup>J</sup>** <sup>=</sup> **<sup>V</sup>***<sup>T</sup> ∏ ∏<sup>T</sup>* **<sup>V</sup>** <sup>=</sup> (**K***H R*)*<sup>T</sup>* **<sup>K</sup>***H R*. If .**K***H R* is full rank, then . **<sup>J</sup>** is symmetric definite positive and.**J**−<sup>1</sup> = (**K***H R*)−<sup>1</sup> (**K***H R*)−*<sup>T</sup>* . Then, both the minimization problem and the hyper-reduced equation have a unique solution. The solution of the minimization problem is:

$$\mathbf{q}^f = \mathbf{V} \mathbf{J}^{-1} \mathbf{V}^T \mathbf{II} \,\mathbf{T}^T \mathbf{q},\tag{2.40}$$

$$\mathbf{V} = \mathbf{V} \left(\mathbf{K}^{HR}\right)^{-1} \mathbf{H}^{T} \mathbf{q},\tag{2.41}$$

$$\begin{aligned} &=\mathbf{V}\left(\mathbf{K}^{RR}\right)^{-1}\mathbf{K}^{RR}\mathbf{q},\\ &=\mathbf{V}\left(\mathbf{K}^{HR}\right)^{-1}\mathbf{V}[\mathcal{P},:]^{T}\mathbf{K}[\mathcal{P},:]\,\mathbf{q},\\ &=\mathbf{V}\left(\mathbf{K}^{HR}\right)^{-1}\mathbf{V}[\mathcal{P},:]^{T}\mathbf{F}[\mathcal{P}],\\ &=\widehat{\mathbf{q}}.\end{aligned} \tag{2.43}$$

$$\mathbf{V} = \mathbf{V} \left(\mathbf{K}^{HR}\right)^{-1} \mathbf{V} [\mathcal{P}, \text{:} \text{:} ]^T \mathbf{F} [\mathcal{P}], \tag{2.43}$$

$$\widehat{\mathbf{q}}.\tag{2.44}$$

As an intermediate result, Eq. (2.42) is the oblique projection.

In linear elasticity, the Céa's lemma holds. Let us denote.**q**◦ <sup>∈</sup> <sup>R</sup>*<sup>N</sup>* the minimizer of the upper bound in Eq. (2.3) related to this lemma:

$$
\begin{aligned}
\text{2. Learning Projection-based Reduced-Order Models.}\\ \|\mathbf{q}^{\circ} = \arg\min\_{\mathbf{q}^{\circ} \in \mathbb{R}^{N}} \|\widetilde{\boldsymbol{u}} - \sum\_{i=1}^{N} \boldsymbol{\varphi}\_{i} q\_{i}^{\star} \|,\\ \|\widetilde{v}^{\circ} = \sum\_{i}^{N} \boldsymbol{\varphi}\_{i} q\_{i}^{\circ}.\end{aligned} \tag{2.45}
$$

$$
\widetilde{w}^{\diamond} = \sum\_{i=1}^{N} \mathfrak{p}\_i q\_i^{\diamond}.\tag{2.46}
$$

The best projection of the minimizer .**q**◦ in the approximation space is denoted by .*γ <sup>P</sup>* :

$$\mathbf{y}^{\mathcal{P}} = \operatorname\*{argmin}\_{\mathbf{g} \in \mathbb{R}^n} (\mathbf{q}^\circ - \mathbf{V}\,\mathbf{g})^T \mathbf{M} (\mathbf{q}^\circ - \mathbf{V}\,\mathbf{g}),\tag{2.47}$$

$$\widetilde{\boldsymbol{v}}^{\mathcal{P}} = \sum\_{i=1}^N \boldsymbol{\varphi}\_i (\mathbf{V}\,\mathbf{y}^{\mathcal{P}})\_i.\tag{2.48}$$

$$
\widetilde{v}^{p} = \sum\_{i=1}^{N} \varphi\_i (\mathbf{V} \ \boldsymbol{\nu}^{p})\_i. \tag{2.48}
$$

Let us introduce an ideal reduced basis .**V**◦ <sup>∈</sup> <sup>R</sup>*N*×*<sup>n</sup>* (It assumes that .*<sup>n</sup>* is an ideal reduced dimension) such that: .**q**◦ = **V**◦ *γ* ◦, and .**V**◦*T***MV**◦ = **I**, where . *Mi j* = ⟨*ϕi*, *ϕ <sup>j</sup>*⟩. Hence.*γ <sup>P</sup>* = **V***T***MV**◦ *γ* ◦.

**Property 2:** In linear elasticity, the upper bound of approximation error is increased by a Chordal distance [101] between. **V** and the ideal reduced basis.**V**◦: .|~*u* −~v◦ |≤ |~*<sup>u</sup>* <sup>−</sup>~v*<sup>P</sup>* |+|*<sup>γ</sup>* ◦

$$\|\widetilde{\boldsymbol{\mu}} - \widetilde{\boldsymbol{\upsilon}}^{\diamond}\| \leq \|\widetilde{\boldsymbol{\mu}} - \widetilde{\boldsymbol{\upsilon}}^{\prime}\| + \|\boldsymbol{\mathcal{y}}^{\diamond}\|\_{2} \, d^{Ch}(\mathbf{V}^{\diamond}, \mathbf{V}),\tag{2.49}$$

where.*dCh*(**V**◦, **V**) is the Chordal distance between.**V**◦ and. **V**.

Hence, the smaller the Chordal distance between the sub-space spanned by. **V** and .**V**◦, the better the reduced prediction by using a Galerkin projection (When the RID covers the full domain). A certification of the reduced projection can be achieved, when all errors admit an upper bound, by following the constitutive relation error proposed in [ 47, 55].

The Chordal distance uses the principal angles .*<sup>θ</sup>* <sup>∈</sup> <sup>R</sup>*<sup>n</sup>*, .θ*<sup>k</sup>* ∈ [0,π/2[ for . *<sup>k</sup>* <sup>=</sup> 1,..., *n*, computed via a full singular value decomposition:

$$\mathbf{V}^{T}\mathbf{M}\mathbf{V}^{\circ} = \mathbf{U}\cos(\theta)\mathbf{U}^{\circ T}, \quad \mathbf{U}^{T}\mathbf{U} = \mathbf{U}^{\circ T}\mathbf{U}^{\circ} = \mathbf{I},\tag{2.50}$$

$$d^{Ch}(\mathbf{V}^{\circ}, \mathbf{V}) = \|\sin(\theta)\|\_{F},\tag{2.51}$$

$$\|\mathbf{U}^{\circ}\|\_{F}^{2} = n,\tag{2.52}$$

where .|·|*<sup>F</sup>* is the Frobenius norm. Here, .cos(*θ*) and .sin(*θ*) are cosine and sine diagonal matrices. In addition the following property holds when a full SVD is computed:

$$\mathbf{U}^{\circ} \mathbf{U}^{\circ T} = \mathbf{I}, \quad \mathbf{U} \, \mathbf{U}^{T} = \mathbf{I}. \tag{2.53}$$

The proof of the previous property is straightforward by using the triangular inequality. We just need to prove that: .|~v◦ <sup>−</sup>~v*<sup>P</sup>* |≤|*<sup>γ</sup>* ◦

$$\|\widetilde{v}^{\diamond} - \widetilde{v}^{P}\| \le \|\mathfrak{p}^{\diamond}\|\_2 \, d^{Ch}(\mathbf{V}^{\diamond}, \mathbf{V}).$$

Hence, the proof is the following:

$$\begin{aligned} \text{the proof is the following:}\\ \|\|\widetilde{v}^{\diamond} - \widetilde{v}^{p}\|\|^{2} &= \mathbf{y}^{\diamond T} \left(\mathbf{V}^{\diamond} - \mathbf{V}\,\mathbf{V}^{T}\mathbf{M}\mathbf{V}^{\diamond}\right)^{T} \mathbf{M} \left(\mathbf{V}^{\diamond} - \mathbf{V}\,\mathbf{V}^{T}\mathbf{M}\mathbf{V}^{\diamond}\right) \mathbf{y}^{\diamond} \\ &= \mathbf{y}^{\diamond T} \left(\mathbf{I} - \mathbf{V}^{\diamond}\mathbf{M}\mathbf{V}\,\mathbf{V}^{T}\mathbf{M}\mathbf{V}^{\diamond}\right) \mathbf{y}^{\diamond}. \end{aligned} \tag{2.54}$$

Therefore

$$
\hat{\mathbf{v}}^{\circ} - \hat{\mathbf{v}}^{\circ} \hat{\mathbf{I}}^{\circ} = \mathbf{y}^{\circ T} \left( \mathbf{I} - \mathbf{U}^{\circ} \cos(\theta) \mathbf{I}^{\circ T} \right) \mathbf{y}^{\circ} \tag{2.55}
$$

$$
\hat{\mathbf{v}}^{\circ} = \mathbf{U}^{\circ} \mathbf{U}^{\circ} - \mathbf{U}^{\circ} \cos(\theta) \mathbf{I}^{\circ T} \mathbf{U}^{\circ T} \tag{2.55}
$$

$$\mathbf{J} = \mathbf{y}^{\circ T} \underbrace{\mathbf{U}^{\circ} (\mathbf{I} - \cos(\theta)^2)}\_{} \mathbf{U}^{\circ T} \mathbf{y}^{\circ} \tag{2.56}$$

$$=\mathbf{y}^{\circ^{T}}\,\mathbf{U}^{\circ}\sin(\theta)^{2}\,\mathbf{U}^{\circ^{T}}\,\mathbf{y}^{\circ}\tag{2.57}$$

$$= \|\sin(\theta)\,\mathbf{U}^{\circ T}\,\mathbf{y}^{\circ}\|\_{2}^{2}.\tag{2.58}$$

For all matrices.**<sup>A</sup>** <sup>∈</sup> <sup>R</sup>*n*×*<sup>m</sup>* and.**<sup>B</sup>** <sup>∈</sup> <sup>R</sup>*m*×*<sup>n</sup>* the following property holds:

$$\|\mathbf{A}\mathbf{B}\|\_{F} \le \|\mathbf{A}\|\_{F} \|\mathbf{B}\|\_{F},$$

and for.**<sup>a</sup>** <sup>∈</sup> <sup>R</sup>*<sup>n</sup>*:.|**a**|*<sup>F</sup>* = |**a**|2. Thus:

$$\begin{array}{cccc} \cdots & \cdots & \cdots & \cdots\\\|\mathbf{R}^{\circ} \colon \|\mathbf{a}\|\_{F} = \|\mathbf{a}\|\_{2} \,. &\\\\|\|\widetilde{v}^{\circ} - \widetilde{v}^{p}\|^{2} \le \|\sin(\theta)\|\_{F}^{2} \, \|\mathbf{U}^{\circ T} \: \mathbf{y}^{\circ}\|\_{2}^{2} \le \|\sin(\theta)\|\_{F}^{2} \, \|\mathbf{y}^{\circ}\|\_{2}^{2} .\end{array} \tag{2.59}$$

**Property 3:** When the identity matrix is substituted for. **K** in Eqs. (2.35) and (2.36) is known as the Gappy POD reconstruction [ 35] of truncated variables .**q**[P]. The reconstructed vector in.R*<sup>N</sup>* is: .~**<sup>q</sup>** <sup>=</sup> **<sup>V</sup>** (**V**[P, :]*<sup>T</sup>* **<sup>V</sup>**[P, :])

$$\widetilde{\mathbf{q}} = \mathbf{V} \left( \mathbf{V}[\mathcal{P}, \colon]^T \mathbf{V}[\mathcal{P}, \colon] \right)^{-1} \mathbf{V}[\mathcal{P}, \colon]^T \mathbf{q}[\mathcal{P}].\tag{2.60}$$

This Gappy POD reconstruction is useless for displacement variables because the oblique projection in Eq. (2.36) is a direct outcome of the hyper-reduced prediction. But such a reconstruction is very convenient for stress variables that the hyperreduced scheme forecasts only on.Ω*R*. The reconstructed stress variables reads: .~**q**<sup>σ</sup> <sup>=</sup> **<sup>V</sup>**<sup>σ</sup> (**V**<sup>σ</sup> [P<sup>σ</sup>

$$\widetilde{\mathbf{q}}^{\sigma} = \mathbf{V}^{\sigma} \left( \mathbf{V}^{\sigma} [\overline{\mathcal{P}}^{\sigma}, :]^{T} \mathbf{V}^{\sigma} [\overline{\mathcal{P}}^{\sigma}, :] \right)^{-1} \mathbf{V}^{\sigma} [\overline{\mathcal{P}}^{\sigma}, :]^{T} \mathbf{q}^{\sigma} [\overline{\mathcal{P}}^{\sigma}], \tag{2.61}$$

where .P<sup>σ</sup> is the set of all stress indices available in .Ω*R*. Since the RID contains interpolation points for.**V**<sup>σ</sup> , these points are included in.P<sup>σ</sup> , therefore the truncated matrix.**V**<sup>σ</sup> [P<sup>σ</sup> , :] is full column rank and the reconstruction is a well posed problem.

**Remark about the RID construction and the DEIM:** If the RID contains solely the elements connected to interpolation points related to the reduced basis . **V**, such that.P = P*<sup>u</sup>*, then the Gappy POD gives the interpolation scheme of the DEIM: .~**q***DEIM* <sup>=</sup> **<sup>V</sup>** (**V**[P*<sup>u</sup>*, :])

$$
\widetilde{\mathbf{q}}^{DEIM} = \mathbf{V} \left( \mathbf{V}[\mathcal{P}^{\mu}, \colon \mathbf{j}]^{-1} \mathbf{q}[\mathcal{P}^{\mu}].\right.\tag{2.62}
$$

But, when considering the hyper-reduction scheme, one can observe overfitting in the sense that the train set of displacement is very well approximated by the DEIM reconstruction, but the hyper-reduced predictions are not accurate. For this reason, we recommend the use of the additional reduced basis.**V**<sup>σ</sup> and the related interpolation points.

Various applications of the hyper-reduction method using a RID have been developed for:


.


### *2.3.5 Hyper-Reduction via Empirical Cubature*

∑

To assemble the linearized equations of the reduced Newton algorithm (2.13) when using the ROM in the online phase, hyper-reduction techniques via empirical cubature aim to compute the costly integrals over the high-fidelity domain by replacing the high-dimensional quadrature formula by a low-dimensional reduced quadrature with positive weights. The ECSW [ 36], ECM [ 45] and LPEQP [100] are methods implementing such reduced quadratures. In this section, we present the ECM, more details are available in [ 19].

We consider the high-fidelity model described in Sect. 2.2. The integrals involved in the assembling of the linearized Eq. (2.7) make use of high-fidelity quadrature formulas. Apply such quadrature to the reduced internal forces vector: *i* (*t*) := { σ ( )

)

$$\begin{split} \hat{F}^{\text{int}}\_{i}(t) &:= \int\_{\Omega} \sigma\left(\epsilon(\hat{u}), \mathbf{y}\right)(\mathbf{x}, t) : \epsilon\left(\psi\_{i}\right)(\mathbf{x}) \\ &= \sum\_{\epsilon \in E} \sum\_{k=1}^{n\_{\epsilon}} \omega\_{k} \sigma\left(\epsilon(\hat{u}), \mathbf{y}\right)(\mathbf{x}\_{k}, t) : \epsilon\left(\psi\_{i}\right)(\mathbf{x}\_{k}), 1 \le i \le n, \end{split} \tag{2.63}$$

where. *E* denotes the set of elements of the mesh,.*ne* the number of quadrature points for the element . *e*, .ω*<sup>k</sup>* and .*xk* are the quadrature weights and points associated to . *e*. The total number of quadrature points is denoted.*NG*.

The ECM aims to approximate the high-fidelity quadrature by a reduced quadrature with positive weights, which, when applied to the reduced internal forces vector, writes

2.3 Linear Manifold Learning for Projection-Based … 25

$$\text{earmanifold Learning for Projection-Bascd}\dots\tag{25}$$

$$\hat{F}\_i^{\text{int}}(t) \approx \sum\_{k'=1}^{n\_\delta} \hat{\alpha}\_{k'} \sigma\left(\epsilon(\hat{u}), \mathbf{y}\right) (\hat{\mathbf{x}}\_{k'}, t) : \epsilon\left(\psi\_i\right)(\hat{\mathbf{x}}\_{k'}), 1 \le i \le n,\tag{2.64}$$

where .ωˆ *<sup>k</sup>*' > 0 and .*x*ˆ*<sup>k</sup>*' are respectively the reduced quadrature weights and points, and.*ng* < *NG* is the length of the reduced quadrature. Denote . *fq* := <sup>σ</sup> (

)

∈(*u*(*q*//*n*)+1), *y* : ∈ ψ(*<sup>q</sup>*%*n*)+<sup>1</sup> , .1 ≤ *q* ≤ *nNc*. where .// and . % are the quotient and the remainder of the Euclidean division. Denote as well . Z*nG* a subset of .[1; *NG*] of size .*nG* and .*J*Z*nG* <sup>∈</sup> <sup>R</sup>*nNc*×*nG* and .*<sup>b</sup>* <sup>∈</sup> <sup>N</sup>*nNc* such that for all .1 ≤ *q* ≤ *nNc* and all.1 ≤ *q*' ≤ *nG*, ()({)

$$J\_{\mathfrak{L}^{\mathfrak{n}\_{G}}} = \left(f\_{q}(\mathbf{x}\_{\mathfrak{L}\_{q}^{\ast \mathcal{G}}})\right)\_{1 \le q \le nN\_{\mathfrak{c}}, \ q' \in \mathfrak{L}^{\mathfrak{n}\_{G}}}, \qquad \mathbf{b} = \left(\int\_{\mathfrak{Q}} f\_{q}\right)\_{1 \le q \le nN\_{\mathfrak{c}}},\tag{2.65}$$

where .Z*nG <sup>q</sup>*' denotes the .*q*' -th element of .Z*nG* . We remind that . *n* is the number of POD modes, see Sect. 2.3.3.2. Let.*ω*<sup>ˆ</sup> <sup>∈</sup> <sup>R</sup>+*<sup>n</sup> G*,. ( *J*Z*nG ω*ˆ *q* = ∑*nG q*' =1 ωˆ *q*'σ ( ∈(*u*(*q*//*n*)+<sup>1</sup>), *y* (*x*Z*nG <sup>q</sup>*' ) : <sup>∈</sup> ( ψ(*<sup>q</sup>*%*n*)+<sup>1</sup> ) (*x*Z*q*'), .1 ≤ *q* ≤ *nNc*, is a candidate approximation for . { Ω σ ( ∈(*u*(*q*//*n*)+<sup>1</sup>), *y* ) : ∈ ( ψ(*<sup>q</sup>*%*n*)+<sup>1</sup> ) = *bq* , .1 ≤ *q* ≤ *nNc*. The problem of finding ()||||

the more accurate reduced quadrature formula of length.*nG* for the reduced internal forces vector is:

$$\left(\hat{\boldsymbol{\omega}}, \mathcal{Z}^{n\_{\mathcal{O}}}\right) = \arg\min\_{\hat{\boldsymbol{\alpha}}' \in \mathbb{R}^{+\pi\_{\mathcal{O}}}, \mathcal{Z}^{n\_{\mathcal{O}}} \subset \left\| 1; N\_{\mathcal{O}} \right\|} \left\| \boldsymbol{J}\_{\mathcal{Z}^{n\_{\mathcal{O}}}} \hat{\boldsymbol{\alpha}}' - \boldsymbol{b} \right\|, \tag{2.66}$$

where.|·| denotes the Euclidean norm. Minimizing the length of the reduced quadrature formula as well leads to a NP-hard problem, which solution can be approximated using a Nonnegative Orthogonal Matching Pursuit algorithm, see Algorithm 2.2.

### **Algorithm 2.2:** Nonnegative Orthogonal Matching Pursuit.

.

**Input :** *J* , *b*, tolerance ∈Op.comp., *xk* ,1 ≤ *k* ≤ *NG*  **Output:** ωˆ *<sup>k</sup>* , *x*ˆ*<sup>k</sup>* ,1 ≤ *k* ≤ *d*. **<sup>1</sup>**Set Z = ∅, *k*' = 0, ωˆ = 0 and *r*0 = *b* ; // initialization **2 while** |*rk*'|2 > ∈Op.comp. |*b*|2 **do <sup>3</sup>** <sup>Z</sup><sup>←</sup> <sup>Z</sup><sup>∪</sup> max index (( *J*[1;*NG* ] )*T rk*' ) **<sup>4</sup>** ωˆ ← arg ωˆ'>0 min | |*b* − *J*Z ωˆ' | | 2 **<sup>5</sup>** *rk*'+1 ← *b* − *J*Z ωˆ **<sup>6</sup>** *k*' ← *k*' + 1 **7 end <sup>8</sup>***d* ← *k*' **<sup>9</sup>***x*ˆ*k* := *x*Z*k* ,1 ≤ *k* ≤ *d*

In Algorithm 2.2,.*J*[1;*NG* ] satisfies the definition (2.65) with.Z*nG* = [1; *NG*]. The positivity of the weights of the reduced quadrature preserves the spectral properties of the operator associated with the high-fidelity problem, see [ 19, Remark 1].

### *2.3.6 Computational Complexity*

In this section we restrict our attention to elliptic problems or to linearized problems. The bilinear part of the weak form for finite-dimensional solution spaces is a matrix. When using Finite Element solution space, this matrix is sparse. But when using a reduced solution space, this matrix is usually a full matrix. Therefore, computational complexity of the finite element prediction is the complexity of the solution of sparse linear system. It scales linearly with.*N* ω2, where . ω is the band width of the sparse matrix. For the reduced prediction, the solution of a full linear system scales linearly with. *n*3. We recommend to restrict linear model reduction schemes, with or without hyper-reduction, to reduced dimension. *n* lower than.*N*1/3, otherwise the solution of reduced equation will have a computational complexity similar to the complexity of the finite element model. This recommendation does not concern explicit solvers.

### **2.4 Nonlinear Manifold Learning for Projection-Based Reduced-Order Modeling**

Consider a parametrized variability, and a set of snapshots generated using the highfidelity model over a sampling of the parameter domain. The parametrized problem is said nonreducible when applying a linear data compression over this set of snapshots leads to a ROB containing too many vectors for the online problem to feature an interesting speedup. Formally, this happens when the Kolmogorov n-width . *dn*(M) decreases too slowly with respect to . *n*, where we recall that . *n* is the cardinality of the ROB,

$$d\_n(\mathcal{M}) := \inf\_{\mathcal{H}\_n \in \text{Gr}(n, \mathcal{H})} \sup\_{u \in \mathcal{M}} \inf\_{v \in \mathcal{H}\_n} ||u - v||\_{\mathcal{H}},\tag{2.67}$$

with the Grassmannian .Gr(*n*, H) being the set of all .*n*-dimensional subspaces of .H and.H*<sup>n</sup>* ∈ Gr(*n*, H) the subspace spanned by the considered ROB. Qualitatively, the solution manifold .M covers too many independent directions to be embedded in a low-dimensional subspace. To address this issue, several techniques have been developed:

• Problem-specific methods tackle the difficulties of some specific physics problems that are known to be nonreducible, such as advection-dominated problems which have been largely investigated, for instance in [ 16, 49, 85].


In the following Sects. 2.4.1 and 2.4.2, we provide more details on the last two entries of the previous list.

### *2.4.1 Nonlinear Dimensionality Reduction via Auto-Encoder*

Nonlinear manifold learning means that the solution manifold is approximated by a domain in the ambient solution space that is not included in a low-dimensional vector subspace, as illustrated in Fig. 2.3.

Let us consider a formal representation of parabolic and nonlinear Partial Differential Equations (PDEs) that are parameterized with respect to some physical parameters of interest. We mean by physical parameters the parameters that appear directly within the equations such as the boundary conditions, the viscosity for fluid mechanics (henceforth the Reynolds number), the time step for dynamical systems of fluid flows or infectious diseases, etc. These parameters are denoted. μ without any loss of generality as introduced in the preceding section. The formal representation of the equations is given as follows: ∂~*u*

.

the preceding section. The formal representation is:

$$
\frac{\partial \widetilde{u}}{\partial t} = f(\widetilde{u}, \mu). \tag{2.68}
$$

**Fig. 2.3** Nonlinear manifold learning

Reduced order modeling based on nonlinear data compression techniques might be a solution for example when the described physical fields by the model equations require a large number of vectors in the ROM as specified above. Nevertheless there are cases where even if the physical solution fields are completely reducible, the Galerkin projection may not be appropriate for the model equations describing this physics.

Convection-diffusion PDEs have this stability issue even more generally when considering Galerkin projection using finite element basis functions. In the literature, it is proved that the coherent structures of a turbulent, unsteady and in-compressible fluid flow are reproducible by a small number of POD basis functions. However, if these functions are used to solve the Newton-Raphson problem of the associated reduced order Galerkin dynamical system, then an instability appears as a function of the time. In the literature, many solutions are proposed to tackle this difficulty while keeping the reduced order approximation in a linear space spanned by the POD basis functions. We can refer to the Petrov-Galerkin technique, the least square minimization of the equations residual, the variational finite element method, etc. Recently, nonlinear approximations of the solution fields in a manifold of reduced dimension start to gain importance in the literature. In this case, the reduced order model is said to be a nonlinear projection based reduced model. Some authors introduced nonlinear approximations using Deep Learning approaches for projection based reduced models. We find in the literature more classical nonlinear approximations based on the Kernel POD technique.

In what follows, we make a focus from the literature on Deep Learning projection based reduced models and their different possible formulations. More precisely, we are talking about Deep AutoEncoders (DAEs) from the domain of Deep Learning. DAEs are artificial neural networks formed of layers of spatial convolutions, nonlinear activation functions and linear systems called fully connected functions. These architectures are used to perform nonlinear dimensional reduction, following unsupervised data compression. Henceforth, the DAE allows to determine latent features within a set of given inputs.

We denote by. *h* and. *g* respectively the encoder mapping and the decoder mapping of a DAE. In general . *g* is the transpose mapping of . *h*. We denote by . α the reduced latent features inferred by . *h*. The dimension of the latent features is equal to the intrinsic dimensionality of the manifold as stated in Remark 2.1 in [ 63]. This intrinsic dimensionality is in the current case the dimension .*N*<sup>μ</sup> of the vector of parameters . μ, which may include the time variable also.

We note that the reduced latent features of a DAE are not parameterized variables in general. In other words they can be seen as non-parametric features associated with a given set of inputs. This formulation is interesting in the framework of projectionbased model reduction, where the associated parameters are given straightforwardly by the physical equations. Hence, knowing the variable parameters within the inputs data helps only with the determination of the intrinsic dimension of the manifold.

In the literature, we find two different formulations for DAEs projection based reduced models.

The first formulation was proposed by Kashima [ 53] and Hartman et al. [ 43]. It is very analogical to the Galerkin formulation of projection based reduced models: given . *<sup>h</sup>* and . *<sup>g</sup>* such that .*<sup>g</sup>* ◦ *<sup>h</sup>* : <sup>~</sup>*<sup>u</sup>* −→ *<sup>u</sup>* and . *<sup>u</sup>* <sup>≈</sup> <sup>~</sup>*u*, then the reduced model is formulated as follows: α(μ) <sup>=</sup> *<sup>h</sup>* ◦ *<sup>f</sup>* ◦ *<sup>g</sup>*( α(μ)) , α(*<sup>t</sup>* <sup>=</sup> <sup>0</sup>) <sup>=</sup> *<sup>h</sup>*(~*u*(*<sup>t</sup>* <sup>=</sup> <sup>0</sup>)) (2.69)

$$\begin{split} \frac{\partial}{\partial t} \widehat{\alpha}(\mu) &= h \circ f \circ \operatorname{g}(\widehat{\alpha}(\mu)) \quad , \quad \widehat{\alpha}(t=0) = h(\widetilde{\mu}(t=0)) \\ &= h \left( f \left( \operatorname{g}(\widehat{\alpha}(\mu)), \mu \right) \right) . \end{split} \tag{2.69}$$

$$\hat{\rho} = h\left(f\left(\mathbf{g}(\widehat{\boldsymbol{\alpha}}(\mu)), \boldsymbol{\mu}\right)\right). \tag{2.70}$$

In the above formulation, the authors relied on the following three points in order to set the time derivative of the latent features equal to the right hand side of Eq. (2.69). *<sup>g</sup>*( α(μ)) belongs to the manifold described by.<sup>μ</sup> −→ <sup>~</sup>*u*(μ), (

• . ∂ ∂*t* • .*h* ∂ ∂*t g*( α(μ))) = ∂ ∂*t <sup>h</sup>* (*g*( α(μ))), • .*h* ◦ *g* = *IN*<sup>μ</sup> .

.

.

**Remark 2.1** The first two above items are hypothesis that are fulfilled in the case where. *h* and. *g* are linear or affine functions. The last item is fulfilled theoretically by the inputs data compression using parameters optimisation of the DAE architecture.

The second formulation of DAEs projection based reduced models is proposed in [ 63], where a least square minimisation of the residual of parabolic PDEs because of the decoder approximation is performed. Then, the reduced model is formulated as follows in order to determine the reduced latent features: α(μ) <sup>=</sup> argmin v(μ)∈R*N*<sup>μ</sup> <sup>|</sup>*<sup>J</sup>* ( *<sup>a</sup>*(μ)) v(μ) <sup>−</sup> *<sup>f</sup>* (*g*( *<sup>a</sup>*(μ)), μ)|<sup>2</sup>

$$\frac{\partial}{\partial t}\widehat{a}(\mu) = \operatorname{argmin}\_{\widehat{v}(\mu)\in\mathbb{R}^{N\mu}} \|J(\widehat{a}(\mu))\widehat{v}(\mu) - f(g(\widehat{a}(\mu)), \mu)\|\_{2}^{2},\tag{2.71}$$
 
$$\text{where } \widehat{v}(\mu) \text{ is the time derivative of } \widehat{a}(\mu), \|.\|\_{2} \text{ denotes the mean square norm or the variance of } \widehat{v}(\mu).$$

euclidean norm and,. *J* is the Jacobian matrix of the decoder mapping which belong

element-wise to the tangent space to the solutions manifold at a given point. . *J* is expressed as follows: . *<sup>J</sup>* : *<sup>a</sup>* −→ ∇*g*( *<sup>a</sup>*).

$$J: \widehat{a} \to \nabla \widehat{g}(\widehat{a}).$$

In this second formulation, the authors do not suppose that the velocity of the decoder approximation is in the manifold of the solutions because mathematically it belongs to the tangent space to the manifold at a given point. Hence, they claim that encoding the decoder approximation will produce a poor approximation by the reduced model.

### *2.4.2 Piecewise Linear Dimensionality Reduction via Dictionary-Based ROM-Nets*

Parts of this section has been inspired from the authors previous work [ 30].

Piecewise linear manifold learning means that the solution manifold is approximated by a dictionary of local linear subspaces, as illustrated in Fig. 2.4, where we denote.M the solution manifold.

The solution manifold is partitioned to get a collection of subsets.M*<sup>k</sup>* ⊂ M that can be covered by a dictionary of low-dimensional subspaces, enabling the use of local linear ROMs. If.{M*<sup>k</sup>* }*k*∈[[1;*K*]] is a partition of.M, then:

$$\forall k \in \llbracket 1; K \rrbracket, \ \forall N \in \mathbb{N}^\*, \quad d\_N(\mathcal{M}\_k) \le d\_N(\mathcal{M}).\tag{2.72}$$

The concepts of ROM-net and dictionary-based ROM-net are introduced in [ 27], which we present in this section. Suppose we dispose of an already computed dictionary of ROMs for the parametrized problem (2.4), where each element of the dictionary is a ROM that can approximate the problem on a subset of the solution manifold .M. A dictionary-based ROM-net is a machine learning algorithm trained to assign the parameter.μ ∈ X to the ROM of the dictionary leading to the most accu-

**Fig. 2.4** Piecewise linear manifold learning

**Fig. 2.5** Exploitation phase of a dictionary-based ROM-net. .*K* local ROMs, combined with a classifier .C*<sup>K</sup>* for automatic ROM .C*<sup>K</sup>* (μ) recommendation, are used to predict the quantity of interest. *Z*(μ)

rate reduced prediction. This assignment, called model recommendation in [ 77], is a classification task, see Fig. 2.5.

The dictionary of ROMs is constructed in a clustering stage, during which snapshots are regrouped depending on their respective proximity on .M, in the sense of a particular dissimilarity measure we introduced in [ 29] and [ 26]. The dissimilarity between two parameter values.μ, μ' ∈ X, denoted by.δ(μ, μ' ), involves the sine of the principal angles between subspaces associated to the solutions of the HFM .*u*(μ), *u*(μ' ) ∈ M, see [ 26, Definition 4.10]. Applying a k-medoids clustering algorithm on the solution manifold.M using the dissimilarity. δ leads to an optimal partitioning for a dictionary of local ROMs, in a sense introduced in [ 26, Property 4.13]. We refer to the remaining of [ 26] for the description of a practical efficiency criterion of the dictionary-based ROM-net, which enables to decide, before the computationally costly steps of the workflow, if a dictionary of ROMs is preferable to one global ROM, and how to calibrate the various hyperparameters of the ROM-net.

**Remark 2.2** Importance of the classification. One could argue that the classification step can be replaced by choosing the cluster . *k* for which the dissimilarity measure .δ(μ,μ˜ *<sup>k</sup>* ) between the parameter. μand the cluster medoid.μ˜ *<sup>k</sup>* is the smallest. However, we recall that the computation of the dissimilarity measure requires solving the HFM at the parameter value . μ, which would render the complete model reduction framework useless. Hence, the classification step enables to bypass this HFM solve and directly recommend the appropriate local ROM.

As briefly mentioned in the introduction of Sect. 2.4, local ROMs can be constructed by partitioning the parameter space [ 33, 34], in which case the classification step is not required: the cluster affectation is made by computing distances directly in the parameter space. In other cases, partitioning in the solution space can be considered without requiring a classification step [ 10]. Consider a time-dependent problems where the initial condition is not a parameter of the problem, and suppose an efficient computation of the clustering distance in the solution space based on the reduced solution at the previous time-step. Then, local basis affectation and switching is possible without requiring classification.

The training of the classifier can be difficult when working with physical fields: simulations are costly, data are in high dimension and classical data augmentation techniques for images cannot be applied. Hence, we can consider replacing the HFM by an intermediate-fidelity solver for generating the data needed for the training of the classifier, by considering coarser meshes and fewer time steps. We point out that the HFM should be used at the end for generating the data required in the training of the local ROMs. We propose in [ 28] improvements for the training of the classifier in our context by developing a fast variant of the mRMR [ 83] feature selection algorithm, and new class-conserving transformations of our data, acting like a data augmentation procedure.

### **2.5 Iterative and Greedy Strategies**

For the sake of the presentation, we have separated the offline and online phases, where the Reduced-Order Model is learnt, then exploited. Actually, more involved strategies exist, where the ROM is constructed in a iterative fashion. The Reduced Basis Method [ 84] is a greedy method, where the ROB is constructed by a single snapshot, corresponding to a randomly chosen parameter value, and the ROB is enriched by the parameter value that maximizes the error made by the current ROM. In complex parameter dependencies, the hyper-reduction scheme can be simulateously constructed as the ROB grows, see [ 31] for such a scheme, with the EIM as hyper-reduction. This greedy construction as been extended to time-dependant problems in the POD-greedy method [ 41], and simultaneous hyper-reduction construction have also been proposed [ 88].

Such iterative strategies rely on a efficient computation of the error made by the ROM. Error estimation is investigated in the next chapter.

### **References**


**Open Access** This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

### **Chapter 3 Error Estimation**

### **3.1 Confidence and Trust in Model-Based Engineering Assisted by AI**

Consider first data-based machine learning techniques. They rely on large sets of examples provided during the training stage and do not learn with equations. Dealing with a situation that do not belong to the training set variability, namely an outof-distribution sample, can be very challenging for these techniques. Trusting them could imply being able to guarantee that the training set covers the operational domain of the system to be trained. Besides, data-based AI can lack in robustness: examples have been given of adversarial attacks in which a classifier was tricked to infer a wrong class only by changing a very small percentage of the pixels of the input image. These models often also lack explainability: it is hard to understand what is exactly learned, what phenomenon occurs through the layers of a neural network. In some cases, information on the background of a picture is used by the network in the prediction of the class of an object, or bias present in the training data will be learned by the AI model, like gender bias in recruitment processes. Addressing these limitations is the subject of the Program Confiance.ai, <sup>1</sup> regrouping French academic as well as industrial partners from defense, transports, manufacturing and energy sectors.

Model-based engineering enjoys better explainability–since the reference-model equations are known and used, and robustness–when the reference-model is wellposed. Concerning trust, in our projection-based reduced-order modeling context, the prediction is in general deterministic, and strict error bounds can be derived in particular cases. This is a main difference with AI-based models, which use notions like confidence intervals and predictive uncertainties, expressed as probability results. In

<sup>1</sup> https://www.confiance.ai/.

<sup>©</sup> The Author(s) 2024

D. Ryckelynck et al., *Manifold Learning*, SpringerBriefs in Computer Science, https://doi.org/10.1007/978-3-031-52764-7\_3

the remainder of this chapter, we present error bounds and indicators in projectionbased reduced-order modeling, depending on the complexity of the underlying equations.

### **3.2 In Linear Elasticity and for Linear Problems**

Parts of this section has been inspired from the authors previous works [ 7] and [ 9].

We suppose that the problem of interest has the following discrete variational form, depending on a parameter. μ in a parameter set. P: for a finite-dimensional space . V of dimension.*N* (with.*N* > 1 resulting, e.g., from discretization), find.*u*<sup>μ</sup> ∈ V such that

$$E\_{\mu}: a\_{\mu}(u\_{\mu}, v) = b(v), \qquad \forall v \in \mathcal{V}, \tag{3.1}$$

where .*a*<sup>μ</sup> is an inf-sup stable bounded sesquilinear form on .V × V and . *b* is a continuous linear form on . V. We define the Riesz isomorphism . *J* from .V' to . V such that for all .*l* ∈ V' and all .*u* ∈ V, .(*Jl*, *u*)<sup>V</sup> = *l*(*u*), where .(·, ·)<sup>V</sup> denotes the inner product of .V with associated norm .|·|<sup>V</sup> and .V' the dual of . V. We denote .βμ := inf *<sup>u</sup>*∈Vsup v∈V |*a*μ(*u*, v)| |*u*|V|v|<sup>V</sup> > 0 the inf-sup constant of.*a*<sup>μ</sup> and.β˜ <sup>μ</sup> a computable positive lower bound of .βμ. For simplicity, we consider that the linear form . *b* is independent of the parameter. μ. The extension to.μ-dependent. *b* is straightforward.

Applying the Galerkin method on the linear problem (3.1), using a ROB .(ψ*i*)<sup>1</sup>≤*i*≤*<sup>n</sup>* <sup>∈</sup> <sup>R</sup>*n*×*<sup>N</sup>* as done in Sect. 2.3.2 leads to finding.*u*ˆ<sup>μ</sup> <sup>∈</sup> <sup>V</sup>*<sup>n</sup>* such that

$$E\_{\mu}: a\_{\mu}(\hat{u}\_{\mu}, u\_{j}) = b(u\_{j}), \qquad \forall j \in \{1, \ldots, n\}. \tag{3.2}$$

The approximate solution on the ROB is written as (2.16):

$$
\hat{u}\_{\mu} = \sum\_{i=1}^{n} \gamma\_i^{\mu} \psi\_i. \tag{3.3}
$$

We assume that the sesquilinear form.*a*<sup>μ</sup> depends on. μ in an affine way, namely there exist . *<sup>d</sup>* functions .α*<sup>k</sup>* (μ) : <sup>P</sup> <sup>→</sup> <sup>C</sup> and .*<sup>d</sup>* .μ-independent sesquilinear forms . *ak* bounded on.V × V such that

$$a\_{\mu}(u,v) = \sum\_{k=1}^{d} a\_k^{\mu} a\_k(u,v), \qquad \forall u, v \in \mathcal{V}, \tag{3.4}$$

which enables the ROM to be online-efficient.

Under the current assumptions, the following error bound holds (see [ 16, Sect. 4.3.2]): for all.μ ∈ P,

### 3.3 In Nonlinear Mechanics of Materials 41

.

.

$$\begin{split} \|\boldsymbol{u}\_{\mu} - \hat{\boldsymbol{u}}\_{\mu}\|\_{\mathcal{V}} \leq \mathcal{E}\_{1}(\mu) &:= \boldsymbol{\beta}\_{\mu}^{-1} \|\boldsymbol{G}\_{\mu}\hat{\boldsymbol{u}}\_{\mu}\|\_{\mathcal{V}}, \\ &= \tilde{\boldsymbol{\beta}}\_{\mu}^{-1} \left\| \boldsymbol{G}\_{00} + \sum\_{i=1}^{\hat{N}} \sum\_{k=1}^{d} \boldsymbol{\alpha}\_{k}^{\mu} \boldsymbol{\gamma}\_{i}^{\mu} \boldsymbol{G}\_{k} \boldsymbol{u}\_{i} \right\|\_{\mathcal{V}}, \end{split} \tag{3.5}$$

with.*G*<sup>μ</sup> the linear map from.V to.V such that. V ϶ *u* |→ *G*μ*u* := *J* ( *a*μ(*u*, ·) − *b* ) ∈ V;.*G*<sup>μ</sup> inheriting the affine dependence of.*a*<sup>μ</sup> on. μ since, for all.*u* ∈ V,

$$G\_{\mu}u = -Jb + \sum\_{k=1}^{d} \alpha\_{k}^{\mu} Ja\_{k}(u, \cdot) = G\_{00} + \sum\_{k=1}^{d} \alpha\_{k}^{\mu} G\_{k} u,\tag{3.6}$$

where.*G*<sup>00</sup> := −*J b* ∈ V and.*Gku* := *J ak* (*u*, ·) ∈ V for all.*k* ∈ {1,..., *d*}.

The error bound (3.5) can rewritten in an equivalent way as

$$\begin{split} \mathcal{E}\_{2}(\boldsymbol{\mu}) &= \tilde{\boldsymbol{\beta}}\_{\mu}^{-1} \left( (\boldsymbol{G}\_{00}, \boldsymbol{G}\_{00})\_{\mathcal{V}} + 2 \text{Re} \sum\_{i=1}^{\hat{N}} \sum\_{k=1}^{d} \boldsymbol{\gamma}\_{i}^{\mu} \boldsymbol{\alpha}\_{k}^{\mu} (\boldsymbol{G}\_{k} \boldsymbol{u}\_{i}, \boldsymbol{G}\_{00})\_{\mathcal{V}} \\ &+ \sum\_{i,j=1}^{\hat{N}} \sum\_{k,l=1}^{d} \boldsymbol{\gamma}\_{i}^{\mu} \boldsymbol{\alpha}\_{k}^{\mu} \boldsymbol{\gamma}\_{j}^{\*\mu} \boldsymbol{\alpha}\_{l}^{\*\mu} (\boldsymbol{G}\_{k} \boldsymbol{u}\_{i}, \boldsymbol{G}\_{l} \boldsymbol{u}\_{j})\_{\mathcal{V}} \right)^{\frac{1}{2}}, \\ &= \tilde{\boldsymbol{\beta}}\_{\mu}^{-1} \left( \boldsymbol{\delta}^{2} + 2 \text{Re} (\boldsymbol{s}^{\prime} \hat{\boldsymbol{x}}\_{\mu}) + \hat{\boldsymbol{x}}\_{\mu}^{\*\mu} \boldsymbol{S} \hat{\boldsymbol{x}}\_{\mu} \right)^{\frac{1}{2}}, \end{split} \tag{3.7}$$

where .<sup>δ</sup> := |*G*00|V, . *<sup>s</sup>* and .*x*ˆ<sup>μ</sup> are vectors in .C*dn* with components . *sI* := (*Gkui*, *<sup>G</sup>*00)<sup>V</sup> and .(*x*ˆμ)*<sup>I</sup>* := <sup>α</sup><sup>μ</sup> *<sup>k</sup>* <sup>γ</sup> <sup>μ</sup> *<sup>i</sup>* , and . *<sup>S</sup>* is a matrix in .C*dn*,*dn* with coefficients . *SI*,*<sup>J</sup>* := (*Gkui*, *Glu <sup>j</sup>*)<sup>V</sup> (with . *I* and . *J* re-indexing respectively .(*k*,*i*) and .(*l*, *j*), for all . 1 ≤ *k*,*l* ≤ *d* and all.1 ≤ *i*, *j* ≤ *n*). The. *t* superscript denotes the transposition. The vector . *s* and the matrix. *S* depend on the ROB.(ψ*i*)<sup>1</sup>≤*i*≤*<sup>n</sup>* but are independent of. μ, hence can be precomputed; while the vector .*x*ˆ<sup>μ</sup> depends on the reduced basis approximation .*u*ˆ<sup>μ</sup> via the coefficients .<sup>γ</sup> <sup>μ</sup> *<sup>i</sup>* . A lower bound .β˜ <sup>μ</sup> of the stability constant of .*a*<sup>μ</sup> is also computed in complexity independent of .*N* (which is possible, for example, by the Successive Constraint Method, see [ 10, 13]).

We would like to stress that.E1(μ) = E2(μ) (in infinite precision arithmetic): the indices. 1 and. 2 are used to denote two different ways to compute the same quantity. In particular,.E1(μ) is not online efficient, while.E2(μ) is.

### **3.3 In Nonlinear Mechanics of Materials**

The remaining of this section is inspired from the authors work [ 8].

We look here for an efficient error indicator in the context of general nonlinearities and nonparametrized variabilities in nonlinear structural mechanics. The generality of the assumptions and the complexity of the model lead us to search for quantities correlated to the error made by the ROM with respect to the HFM, instead of rigorous error bounds as considered in the previous section.

The problem of interest is the same as described in Sect. 2.2, and we suppose that we have constructed a ROM following the methods described in Sects. 2.3 and 2.3.5, namely used POD for data compression and ECM for operator compression.

The quantity of interest is the accumulated plastic strain over the complete structure. Since this is a dual quantity, the ROM do not provide directly a prediction over the structure, but only at the reduced quadrature points selected by ECM, see Sect. 2.3.5. The Gappy-POD can be used for to recover the dual quantity of interest over the rest of the structure, see [ 8, Algorithms 3 and 4] for a presentation on the present context, and [ 11] for in seminal paper on Gappy-POD.

A quantification for the prediction relative error of the accumulated plastic strain is defined as

$$E^{P}\_{\mu} := \begin{cases} \frac{\|p\_{\mu} - \hat{p}\_{\mu}\|\_{L^{2}(\Omega)}}{\|p\_{\mu}\|\_{L^{2}(\Omega)}} & \text{if} \|p\_{\mu}\|\_{L^{2}(\Omega)} \neq \mathbf{0},\\ \frac{\|p\_{\mu} - \hat{p}\_{\mu}\|\_{L^{2}(\Omega)}}{\max\limits\_{\mu \in \mathcal{P}\_{\text{off}}} \|p\_{\mu}\|\_{L^{2}(\Omega)}} & \text{otherwise}, \end{cases} \tag{3.8}$$

where .*p*<sup>μ</sup> and .*p*˜<sup>μ</sup> are respectively the high-fidelity and reduced predictions for the accumulated plasticity field at the variability . μ, and .Poff. is the set of variabilities encountered during the offline stage. We underline the fact that .*p*˜<sup>μ</sup> is the reduced prediction over the complete structure, hence after applying the Gappy-POD reconstruction.

Define the ROM-Gappy-POD residual as

$$\mathcal{E}\_{\mu}^{p} := \begin{cases} \frac{\|\hat{\boldsymbol{p}}\_{\mu} - \hat{\boldsymbol{p}}\_{\mu}\|\_{2}}{\|\hat{\boldsymbol{p}}\_{\mu}\|\_{2}} & \text{if} \|\hat{\boldsymbol{p}}\_{\mu}\|\_{2} \neq 0, \\\frac{\|\hat{\boldsymbol{p}}\_{\mu} - \hat{\boldsymbol{p}}\_{\mu}\|\_{2}}{\max\_{\mu \in \mathcal{P}\_{\text{off}}} \|\hat{\boldsymbol{p}}\_{\mu}\|\_{2}} & \text{otherwise}, \end{cases} \tag{3.9}$$

where . *p*˜ <sup>μ</sup> is the reduced prediction (after applying the Gappy-POD) taken at the reduced quadrature points (.*p*˜μ,*<sup>k</sup>* = ˜*p*μ(*x*ˆ*<sup>k</sup>* ), .1 ≤ *k* ≤ *m<sup>p</sup>*), . *p*ˆ <sup>μ</sup> is the vector of the accumulated plastic strain as computed by the constitutive law solver at the reduced quadrature points during the online stage, and .|·|<sup>2</sup> denotes the Euclidean norm. Notice that in the general case, . *p*˜ <sup>μ</sup> /= *p*ˆ <sup>μ</sup>: this discrepancy is at the base of our proposed error indicator.

Notice that the relative error.*E <sup>p</sup>* <sup>μ</sup> involves fields and.*L*2-norms whereas the ROM-Gappy-POD residual .E*<sup>p</sup>* <sup>μ</sup> involves vectors of dual quantities in the set of reduced integration points and Euclidean norms. In (3.9),.| *p*˜ <sup>μ</sup> − *p*ˆ <sup>μ</sup>|<sup>2</sup> is the error between the online evaluation of the accumulated plastic strain by the behavior law solver:. *p*ˆ <sup>μ</sup>, and the reconstructed prediction at the reduced integration points. *x*ˆ*<sup>k</sup>* :. *p*˜ <sup>μ</sup>,.1 ≤ *k* ≤ *m<sup>p</sup>*. It is explained in [ 8, Sect. 4.1] that.| *p*˜ <sup>μ</sup> − *p*ˆ <sup>μ</sup>|<sup>2</sup> is also the residual of the least-square optimization involved in the online stage of the Gappy-POD.

Let .*<sup>B</sup>* <sup>∈</sup> <sup>R</sup>*mp*×*<sup>n</sup> <sup>p</sup>* such that .*Bk*,*<sup>i</sup>* <sup>=</sup> <sup>ψ</sup> *<sup>p</sup> <sup>i</sup>* (*x*ˆ*<sup>k</sup>* ), .1 ≤ *k* ≤ *m<sup>p</sup>*, .1 ≤ *i* ≤ *n <sup>p</sup>*, . *K* := {*p*μ, for all possible variabilities μ} and .*d*(*K*, *W*)*L*2(Ω) := sup v∈*K* inf w∈*W* |v − w|*L*2(Ω), with .*W* a finite-dimensional subspace of.*L*2(Ω). The following propositions and corollary are proven in [ 8, Sect. 4.1].

**Proposition 3.1** *There exist two positive constants*.*C*<sup>1</sup> *and*.*C*<sup>2</sup> *independent of*. μ *(but dependent on*.*n p) such that* 

$$\left\|\left\|p\_{\mu}-\tilde{p}\_{\mu}\right\|\_{L^{2}(\Omega)}^{2}\leq C\_{\text{l}}\left\|Bz\_{\mu}-\hat{p}\_{\mu}\right\|\_{2}^{2}+C\_{\text{l}}\left\|p\_{\mu}-\hat{p}\_{\mu}\right\|\_{2}^{2}+C\_{2}d(K,\text{Span}\{\psi\_{i}^{P}\}\_{1\leq i\leq n^{p}})\_{L^{2}(\Omega)}^{2}.\tag{3.10}$$

**Proposition 3.2** *There exist two positive constants* .*K*<sup>1</sup> *and* .*K*<sup>2</sup> *independent of* . μ *such that* 

$$\|\|\tilde{\mathbf{p}}\_{\mu} - \hat{\mathbf{p}}\_{\mu}\|\_{2}^{2} \leq K\_{1} \left\|p\_{\mu} - \tilde{p}\_{\mu}\right\|\_{L^{2}(\Omega)}^{2} + K\_{2} \|\mathbf{p}\_{\mu} - \hat{\mathbf{p}}\_{\mu}\|\_{2}^{2}.\tag{3.11}$$

**Corollary 3.3.1** *Suppose that the reduced solution is exact up to the considered time step at the online variability* . μ*:* .*p*<sup>μ</sup> = ˜*p*<sup>μ</sup> *in* .*L*<sup>2</sup>(Ω)*. In particular, the behavior law solver has been evaluated with the exact strain tensor and state variables at the integration points* . *xk , leading to* .*p*ˆμ(*x*ˆ*<sup>k</sup>* ) = *p*μ(*x*ˆ*<sup>k</sup>* )*,* .1 ≤ *k* ≤ *m<sup>d</sup> . From Proposition 3.2,*  .<sup>|</sup> *<sup>p</sup>*˜ <sup>μ</sup> <sup>−</sup> *<sup>p</sup>*<sup>ˆ</sup> <sup>μ</sup>|<sup>2</sup> <sup>=</sup> <sup>0</sup>*, and*.E*<sup>p</sup>* <sup>μ</sup> = 0*.* 

We observe that in practice, the evaluations of the ROM-Gappy-POD residual .E*p* <sup>μ</sup> (3.9) and the error.*E <sup>p</sup>* <sup>μ</sup> (3.8) are very correlated in our numerical simulations. The idea is to exploit this correlation by training a Gaussian process regressor for the function.E*<sup>p</sup>* <sup>μ</sup> |→ *<sup>E</sup> <sup>p</sup>* <sup>μ</sup>. At the end of the offline stage, we propose to compute reduced predictions at variability values .{μ*i*}1≤*i*≤*Nc* encountered during the data generation step, and the corresponding couples . ( *E p* <sup>μ</sup>*i*, <sup>E</sup>*<sup>p</sup>* μ*i* ) , .1 ≤ *i* ≤ *Nc*. A Gaussian process regressor is trained on these values and we define an approximation function

$$
\mathcal{E}\_{\mu}^{p} \mapsto \mathbf{G} \mathbf{p} \mathbf{r}^{p} (\mathcal{E}\_{\mu}^{p}), \tag{3.12}
$$

for the error.*E <sup>p</sup>* <sup>μ</sup> at variability. μ as the mean plus 3 times the standard deviation of the predictive distribution at the query point.E*<sup>p</sup>* <sup>μ</sup>: this is our proposed error indicator. If the dispersion around the learning data is small for certain values.E*<sup>p</sup>* <sup>μ</sup>, then adding 3 times the standard deviation will not change very much the prediction, whereas for values with large dispersion of the learning data, this correction aims to provide an error indicator larger than the error. We use the GaussianProcessRegressor python class from scikit-learn [ 17]. Notice that although some operations in computational complexity dependent on.*N* are carried-out, we are still in the offline stage, and they are much faster than the resolutions of the large size systems of nonlinear equations (2.8). If the offline stage is correctly carried-out and since.E*<sup>p</sup>* <sup>μ</sup> is highly correlated with the error, only small values for.E*<sup>p</sup>* <sup>μ</sup> are expected to be computed. Hence, in order to train the Gaussian process regressor correctly for larger values of the error, the reduced Newton algorithm (2.13) is solved with a large tolerance.∈ROM Newton = 0.1. We call these

**Fig. 3.1** Workflow for the offline stage with error indicator calibration [ 8]

operations "calibration of the error indication", see Algorithm 3 for a description and Fig. 3.1 for a presentation of the workflow featuring this calibration step.

**Algorithm 3:** Calibration of the error indicator. **Input** : Outputs of the data generation, data compression and operator compression steps of Section 2.3.3 **Output**: Approximation function E*<sup>p</sup>* <sup>μ</sup> |→ Gpr *<sup>p</sup>*(E*<sup>p</sup>* μ) of the error **<sup>1</sup>**Set Z = ∅, *k*' = 0, ωˆ = 0 and *r*0 = *b* ; // initialization **2 for** *i* = 1 ..., *Nc* **do 3** Construct and solve the reduced problem (2.13) with ∈ROM Newton = 0.1 **<sup>4</sup>** Compute the reconstructed plasticity *p*˜<sup>μ</sup>*i* using Gappy-POD and E*<sup>p</sup>* μ*i*  **<sup>5</sup>** Compute the error *E <sup>p</sup>* <sup>μ</sup>*i* using (3.8) **<sup>6</sup>**X ← X ∪ ( E*p*  <sup>μ</sup>*i* , *E <sup>p</sup>* μ*i* ) **<sup>7</sup>**Construct an approximation function E*<sup>p</sup>* <sup>μ</sup> |→ Gpr *<sup>p</sup>*(E*<sup>p</sup>* μ) of the error *E <sup>p</sup>* <sup>μ</sup> using a Gaussian process regression and the data from X

We recall that in model order reduction, the original hypothesis is the existence of a low-dimensional vector space where an acceptable approximation of the highfidelity solution lies. The hypothesis is formalized under a rate of decrease for the Kolmogorov .*n*-width with respect to the dimension of this vector space. The same hypothesis is made when using the Gappy-POD to reconstruct the dual quantities, which are expressed as a linear combination of constructed modes. For both the primal and dual quantities, the modes are computed by searching some low-rank structure of the high-fidelity data. The coefficients of the linear combination for reconstructing the primal quantities are given by the solution of the reduced Newton algorithm (2.13). After convergence, the residual is small, even in cases where the reduced order model exhibits large errors with respect to the high-fidelity reference: this residual gives no information on the distance between the reduced solution and the high-fidelty finite element space.

However, in the online phase of the ROM-Gappy-POD reconstruction (see [ 8, Algorithm 4]), the coefficients.*p*ˆμ,*<sup>k</sup>* (the accumulated plastic strain computed by the constitutive law solver during the online stage) contain information from the highfidelity behavior law solver. Moreover, an overdetermined least-square is solved, which can provide a nonzero residual that implicitly contains this information from the high-fidelity behavior law solver: namely the distance between the prediction from the behavior law and the vector space spanned by the Gappy-POD modes (restricted to the reduced integration points): this is the term.|*Bz*<sup>μ</sup> − *p*ˆ <sup>μ</sup>|<sup>2</sup> in (3.10). Hence, the ability of the online variability to be expressed on the Gappy-POD modes is monitored through the behavior law solver on the reduced integration points. When the ROM is solved for an online variability not included in the offline variabilities, then the new physical solution cannot be correctly interpolated using the POD and Gappy-POD modes: hence, the ROM-Gappy-residual becomes large.

From Proposition 3.2, if .|*Bz*<sup>μ</sup> − *p*ˆ <sup>μ</sup>|<sup>2</sup> = | *p*˜ <sup>μ</sup> − *p*ˆ <sup>μ</sup>|<sup>2</sup> is large, then the global error . | <sup>|</sup>*p*<sup>μ</sup> − ˜*p*<sup>μ</sup> | | *<sup>L</sup>*2(Ω) and/or the error at the reduced integration points .*x*ˆ*<sup>k</sup>* is large, which makes .|*Bz*<sup>μ</sup> − *p*ˆ <sup>μ</sup>|<sup>2</sup> a good candidate for a enrichement criterion for the ROM. A limitation of the error indicator can occur if the online variability activates strong nonlinearities on areas containing no point from the reduced integration scheme, namely through the term.*C*2*d*(*K*, Span{<sup>ψ</sup> *<sup>p</sup> <sup>i</sup>* }1≤*i*≤*<sup>n</sup> <sup>p</sup>* )<sup>2</sup> *<sup>L</sup>*2(Ω) in (3.10).

We recall that the error indicator (3.12) is a regression of the function.E*<sup>p</sup>* <sup>μ</sup> |→ *<sup>E</sup> <sup>p</sup>* μ. In the online phase, we only need to evaluate .E*<sup>p</sup>* <sup>μ</sup> and do not require any estimation for the other terms and constants appearing in Propositions 3.1 and 3.2.

Equipped with an efficient error indicator, we are now able to assess the quality of the approximation made by the reduced order model in the online phase. If the error indicator is too large, an enrichment step occurs: the high-fidelity model is used to compute a new high-fidelity snapshot, which is used to update the POD and Gappy-POD basis, as well as the reduced integration schemes. Notice that for the enrichment steps to be computed, the displacement field and all the state variables of the previous time step need to be reconstructed on the complete mesh.Ω to provide the high-fidelity solver with the correct material state. The workflow for the *online*  stage with enrichment is presented in Fig. 3.2.

We refer to [ 8] for more details on this subject, and detailed numerical applications in nonlinear structural mechanics for this error indicator and its ability to enrich a ROM in the online stage.

Notice that another (noncertified) indicator in nonlinear solid mechanics with internal variables has been proposed in [ 1], aiming to approximate the dual norm of residuals in the same fashion as in the linear case described in Sect. 3.2. For such nonlinear case, rigorous error bounds are not obtained: a gappy-POD-based approximate representation of the stress tensor is used, and the inf-sup constant evaluation has been replaced by a normalization of the residual using the norm of the Riesz elements for the external loading.

**Fig. 3.2** Workflow for the online stage with enrichment [ 8]

### **3.4 In Computational Fluid Dynamics**

In the following section, we present a priori error estimates due to the POD-Galerkin approximation applied to fluid dynamics equations in particular and parabolic nonlinear PDEs in general. It is a theoretical result on the convergence of the POD-Galerkin reduced order model towards the high-fidelity semi-discretized equations in the sense of the spatial variable. The solution of these semi-discretized equations is denoted .~*u<sup>h</sup>* over a time interval .[0, *<sup>T</sup>* ] such that .*<sup>h</sup>* <sup>=</sup> <sup>1</sup> *<sup>M</sup>* and .*M* is the cardinal of a Hilbert basis that is capable of generating the high-fidelity semi-discretized solution at some specified.*M* time instants. This orthonormal basis can be obtained in an a posteriori fashion thanks to the snapshots POD method. It is denoted by.(ψ*i*)*<sup>i</sup>*=1,...,*<sup>M</sup>* .

The following convergence result includes furthermore a discussion on the stability of the Galerkin projection technique applied to parabolic PDEs. This result has been developed and published in the following papers: [ 2– 5].

In the literature, we can find many works around the problem of defining convergent and a priori upper bounds for the POD-Galerkin reduced models for parabolic equations. It is a subject of great interest, if we can quantify efficiently the error of a reduced solution . *u* obtained with an approximation technique of the corresponding high fidelity semi-discretized solution . ~*u<sup>h</sup>*. This problem can be seen as a theoretical confidence interval around a training point of an approximation model (details on some convergence results of the literature should be added). Let us denote by .<sup>Ω</sup> the open and bounded domain of the spatial variable, such that .<sup>Ω</sup> <sup>∈</sup> <sup>R</sup>*<sup>d</sup>* , where .*d* = 1 *or* 2. Let us consider the parabolic PDEs for which the weak formulation in the space of the solution .*u<sup>h</sup>* spanned by the POD basis of cardinal . *M*, is described as follows:


.


The following result is proved:

**Theorem 3.1** *If*.*n* << *M is the dimension of the truncated POD-Galerkin reduced solution* . *u that represent the training point (solution)* .*uh, then we can derive the following a priori upper bound of the* .[*L*<sup>2</sup>(Ω)] *<sup>d</sup>*<sup>−</sup> *error between* . *<sup>u</sup>*(*t*) *and* .~*u<sup>h</sup>*(*t*)*,* . <sup>∀</sup> .*t* ∈ [0, *T* ]*:* 

$$\left\| \left( \widehat{\boldsymbol{u}}^{h} - \widehat{\boldsymbol{u}} \right) (\boldsymbol{t}) \right\|\_{\left[L^{2}(\Omega)\right]^{d}}^{2} \leq f\_{1}(N). \tag{3.13}$$

*Where*. *f*1(*n*) *is the remainder of the sum of a convergent series;*. *f*1(*n*) *converges to*. 0 *when*.*N converges to*. *M. More precisely,*. *f*1(*n*) *is a function of the remainder of the sum of the POD eigenvalues*.(λ*<sup>k</sup>* )*<sup>k</sup>*=1,...,*<sup>M</sup> obtained from the Snapshots POD applied to*.*M temporal snapshots of*.*uh: it is expressed into two different fashions: Either* 

$$f\_1(n) = \left(1 + 2C\_b K + \frac{C\_a^2}{\epsilon}\right) T \sum\_{k=n+1}^{M} \lambda\_k,\tag{3.14}$$

*if*.(∈ − 2β + 6*CbK*) ≤ 0*, for a strictly positive real number*. ∈*, or* 

$$f\_1(n) = \left(T + \left(2C\_b K + \frac{C\_a^2}{\epsilon}\right)T\exp(T[2C\_a K + K^2(1 + C\_a^2)])\right)\sum\_{k=n+1}^{M} \lambda\_k,\tag{3.15}$$

*if*.(∈ − 2β + 6*CbK*) > 0*, for all strictly positive real numbers*. ∈.

The mathematical proof of Theorem 3.1 is based on the properties of the forms . *a* and . *b*, the application of the Young inequality and the Gronwall lemma. For the details of the proof, refer to [ 4].

**Remark 3.1** The result of Theorem 3.1 is applicable in particular for the.1*D* Burgers equation and the .2*D* unsteady and incompressible Navier-Stokes equations, by remarking the following:

• When .*d* = 1,.*H*<sup>1</sup>(Ω) ⊂ *L*∞(Ω): there exists .*C*<sup>1</sup> <sup>∞</sup> <sup>∈</sup> <sup>R</sup>+∗ such that .∀. <sup>v</sup> <sup>∈</sup> *<sup>H</sup>*<sup>1</sup>(Ω), .|v|*<sup>L</sup>*∞(Ω) ≤ *C*<sup>1</sup> <sup>∞</sup> |∇v|*L*2(Ω).


$$\|\boldsymbol{\upsilon}\|\_{L^{4}(\Omega)} \leq C \left\|\boldsymbol{\upsilon}\right\|\_{L^{2}(\Omega)}^{1/2} \left\|\nabla \boldsymbol{\upsilon}\right\|\_{L^{2}(\Omega)}^{1/2} \dots$$

• If we denote by .*<sup>S</sup>* the square matrix of dimension .*<sup>M</sup>* defined by: . *Si*,*<sup>j</sup>* <sup>=</sup> ⟨ ψ*i*, ψ*<sup>j</sup>* ⟩ [*H*1(Ω)]*<sup>d</sup>* , then.∀. <sup>v</sup> in.V*<sup>h</sup>* the solutions space spanned by the complete POD basis. <sup>ψ</sup> we have the following inequality:.|v|[*H*1(Ω)]*<sup>d</sup>* <sup>≤</sup> √|*<sup>S</sup>*| |v|[*L*2(Ω)]*<sup>d</sup>* . Where .|*S*| is the spectral norm of the matrix.*S*. For details, refer to [ 14].

**Remark 3.2** A result of stability with respect to time of the POD-Galerkin reduced model for nonlinear and dissipated parabolic PDEs can be derived from result 3.1. If .μ denotes only the viscosity constant in this particular case (without any loss of generality), then the POD-Galerkin reduced model is stable with respect to time when. μ satisfies the following inequality:

$$
\mu \ge \frac{6C\_b K + \epsilon}{2\beta},
$$

where. ∈ is a strictly positive real number.

.

**Remark 3.3** More particularly, in the case of linear and dissipated parabolic PDEs, the stability condition of the POD-Galerkin reduced model becomes equivalent to .μ ≥ ∈ 2*cp* . The error of the POD-Galerkin reduced model with respect to the high fidelity training point can be estimated exactly as in the Céa lemma applied for elliptic and sesquilinear PDEs: if.μ ≥ 1 then,

$$\left\|(u^h - \widehat{u})(t)\right\|\_{[L^2(\Omega)]^d}^2 \le \left(1 + \frac{C\_a^2}{2\beta}\right)T \sum\_{n=N+1}^M \lambda\_n.$$

Based on the same methodology, we propose in what follows an a priori error estimate and a convergence result for POD-Galerkin reduced model parameter wise. In other words, we show that when the parameters change with respect to the training ones, a confidence interval is obtained around the new test solution. The width of this confidence interval converges to the truncation error of the ROM at the training parameters. This parametric convergence result is formulated as follows:

**Theorem 3.2** *Let us denote by* .*u<sup>h</sup>* <sup>μ</sup><sup>0</sup> *a training solution associated with a training parameter* .μ0*. So we suppose we have only one training point for the POD-Galerkin reduced model, without any loss of generality. If* .ψ<sup>μ</sup><sup>0</sup> *denotes the Hilbert basis obtained from the POD applied to the high fidelity training solution and*. *u*μ,μ<sup>0</sup> *denotes the truncated POD-Galerkin reduced solution that approaches the test point (solution)*.*u<sup>h</sup>* <sup>μ</sup> *in the reduced POD space of dimension*. *n spanned by*.ψ<sup>μ</sup><sup>0</sup> *, then:*

.

$$\left\| \left( \mu\_{\mu}^{h} - \widehat{u}\_{\mu,\mu\_{0}} \right)(t) \right\|^{2} \leq f\_{1}^{\mu\_{0}}(n) \left( 1 + \left\| \mu - \mu\_{0} \right\|^{2a} \right) + f\_{2}^{\mu\_{0}}(n) \left\| \mu - \mu\_{0} \right\|^{a}, \tag{3.16} \right)$$

*where*. *f* μ0 <sup>2</sup> (*n*) *is the remainder of the sum of a convergent series;*. *f* μ0 <sup>2</sup> (*n*) *converges to*. 0 *when*. *n converges to*. *M. More precisely*. *f* μ0 <sup>2</sup> (*n*) *is a function of the remainder of the sum of the orthogonal projection coefficients of the characteristic function* . 1Ω*<sup>x</sup> such that*.Ω*<sup>x</sup> tends to*.Ω *when*. *x tends to*.∂Ω*: it is expressed as follows:* 

$$f\_2^{\mu\_0}(n) = \left( \left\|(u\_{\mu}^h - u\_{\mu\_0}^h)(0)\right\|^{-4} + \frac{B}{A} (\exp(-2At) - 1) \right)^{-1/2} C\_a^2 \left\|\nabla u\_{\mu\_0}^h\right\|^2 \sum\_{k=n+1}^M \langle 1\_{\Omega\_k}, \psi\_k^{\mu\_0} \rangle^2,\tag{3.17}$$

*where* .*A and*.*B are strictly positive constants of which the detailed expressions are given in [ 2].* 

The proof of the above theorem is published in [ 2]. It is based on the properties of the two forms. *a* and. *b*, the application of the Young inequality and the resolution of a nonlinear ordinary differential inequality of Ricatti type.

### **3.5 A Note on Accuracy of a Posteriori Error Bounds and Round-Off Errors**

In this section, we explain why the online-efficient error bound (3.7) may be sensitive to round-off errors.

In computers, real numbers are represented by a finite number of bits, called floating-point representation. Current architectures are optimized for the IEEE 754 double-precision binary floating-point format. Let . *x* and . *y* be real numbers. When computing the operation.*x* + *y*, the result returned by the computer can be different from its theoretical value. Whenever the difference is substantial, a loss of significance occurs. A well-known case of loss of significance is when. *x* and. *y* are almost opposite numbers. Suppose that.*x* = −*y*. We denote by.maxfl(*x* + *y*) the result that the computer returns when the maximal accumulation of round-off errors occurs when computing the summation. There holds

$$|\text{maxfl}(\mathbf{x} + \mathbf{y})| \approx 2\epsilon |\mathbf{x}|,\tag{3.18}$$

where. ∈ is called the machine precision. In double precision,.∈ = 5 × 10−<sup>16</sup> (see [ 12, Sect. 1.2]).

When implementing an algorithm, one should ensure that each step is free of such a loss of significance. In some cases, simply changing the order of the operations can prevent these situations. As an illustration, consider .*x* = 1, .*y* = 1 + 10−7, and the operation .*x* <sup>2</sup> − 2*x y* + *y*2. This is a sum of terms where the first intermediate result in the sum is.14 orders larger than the result. Therefore, a loss of significance is expected. The relative error of this computation is about .8 × 10−4. Computing .(*x* − *y*)2, which is the factorization of the considered operation, leads to a relative error of about .10−9. Thus, the terms of the sum are only . 7 orders larger than the results, leading to a less catastrophic loss of significance. In this specific case, the remedy consists in carrying out the sum before the multiplication. In our projectionbased ROM context, the evaluation of the formula .E<sup>2</sup> suffers from such a loss of significance, as we now explain.

We investigate the influence of round-off errors when computing the error bounds (3.5) and (3.7) for respectively.E1(μ) and.E2(μ). As observed in the previous paragraph, the computation of a polynomial using a factorized form is more accurate than using the developed form, in particular at points close to its roots. Here,. ( β˜ μE2(μ))<sup>2</sup> is a multivariate polynomial of degree. 2 in.*x*ˆ<sup>μ</sup> computed in a developed form, whereas the scalar product .(*G*μ*u*μ, *G*μ*u*μ)<sup>V</sup> used in the computation of .E1(μ) is not developed. The following holds (see [ 9, Proposition 2.2.1] for the proof)

**Proposition 3.3** *Let*.μ ∈ P *and let*.maxfl(β˜ <sup>μ</sup>E*<sup>k</sup>* (μ))*,*.*k* = 1, 2*, denote the evaluation of*.β˜ <sup>μ</sup>E*<sup>k</sup>* (μ) *when the maximum accumulation of round-off errors occurs. There holds* 

$$\begin{aligned} \text{maxfl}(\bar{\beta}\_{\mu}\mathcal{E}\_{1}(\mu)) &\geq 2\delta\epsilon, \\\text{maxfl}(\tilde{\beta}\_{\mu}\mathcal{E}\_{2}(\mu)) &\geq 2\delta\sqrt{\epsilon}, \end{aligned} \tag{3.19}$$

*where*.δ = |*G*00|<sup>V</sup> *and*. ∈ *is the machine precision.* 

.

From this proposition, we notice that the online-efficient formula .E2(μ) suffers from an important loss of significance.

We present below an error estimator proposed in [ 9] that enjoys both accuracy and online-efficiency. Let.σ := 1 + 2*dn* + (*dn*)2. For a given.μ ∈ Ptrial and the resulting .*u*ˆ<sup>μ</sup> <sup>∈</sup> Span{ψ1, ..., ψ*n*} solving the reduced problem (3.2), we define.*X*ˆ(μ) <sup>∈</sup> <sup>C</sup><sup>σ</sup> as the vector with components.(1, *x*ˆ<sup>μ</sup>*<sup>I</sup>* , *x*ˆ<sup>∗</sup> μ*I* , *x*ˆ<sup>∗</sup> μ*I <sup>x</sup>*ˆ<sup>μ</sup>*<sup>J</sup>* ), where.*x*ˆ<sup>μ</sup>*<sup>I</sup>* <sup>=</sup> <sup>α</sup><sup>μ</sup> *<sup>k</sup>* <sup>γ</sup> <sup>μ</sup> *<sup>i</sup>* (we recall that .γ <sup>μ</sup> *<sup>i</sup>* are the coefficients of the reduced solution in the reduced basis, see (3.3), and .αμ *<sup>k</sup>* the coefficients of the affine decomposition of .*a*<sup>μ</sup> in (3.4)), with . 1 ≤ *I*, *J* ≤ *dn* (with.*I* = *i* + *n*(*k* − 1) such that.1 ≤ *i* ≤ *n*,.1 ≤ *k* ≤ *d*, and with. *J* = *j* + *n*(*l* − 1) such that.1 ≤ *j* ≤ *n*,.1 ≤ *l* ≤ *d*). We can write the right-hand side of (3.7) as a linear form in.*X*ˆ(μ) as follows:

$$\delta^2 + 2\text{Re}(s^\prime \hat{x}\_\mu) + \hat{x}\_\mu^{\ast t} S \hat{x}\_\mu = \sum\_{p=1}^\sigma t\_p \hat{X}\_p(\mu),\tag{3.20}$$

where.*tp* is independent of. μ (as. δ,. *s*, and. *S* are independent of. μ) and.*X*ˆ *<sup>p</sup>*(μ) is the .*p*-th component of.*X*ˆ(μ).

Consider the function of two variables .(*p*, μ) |→ *X*ˆ *<sup>p</sup>*(μ), for all . *p* ∈ {1, ..., σ} and all.μ ∈ P. We look for an approximation of this function in the form

$$\forall \mu \in \mathcal{P}, \forall p \in \{1, \ldots, \sigma\}, \,\,\hat{X}\_p(\mu) \approx \sum\_{r=1}^{\hat{\sigma}} \lambda\_r^{\hat{\sigma}}(\mu) \hat{X}\_p(\mu\_r),\tag{3.21}$$

for a certain parameter.σˆ ≤ σ. The Empirical Interpolation Method (EIM) provides a numerical procedure to construct this approximation and to choose the interpolation points (see [ 6, 15]), which leads to the following formula for computing the error bound

$$\mathcal{E}\_3(\mu) := \tilde{\beta}\_{\mu}^{-1} \left( \sum\_{r=1}^{\hat{\sigma}} \lambda\_r^{\hat{\sigma}}(\mu) V\_r \right)^{\frac{1}{2}},\tag{3.22}$$

where .*Vr* = | <sup>|</sup>*G*<sup>μ</sup>*<sup>r</sup> <sup>u</sup>*ˆ<sup>μ</sup>*<sup>r</sup>* | | 2 <sup>V</sup>, and where .λ<sup>σ</sup>ˆ(μ) and .μ*<sup>r</sup>* are provided by EIM, see [ 9, Sect. 3.2] for all the details of this derivation. There holds (see [ 9, Proposition 3.2.1]):

**Proposition 3.4** *The computation of the formula*.E<sup>3</sup> *is well defined, and this formula is online-efficient.* 

Besides,.E<sup>3</sup> involving a linear combination of accurately computed scalar products (see (3.5)), it is not subject to the loss of significance encountered in.E2.

We refer to [ 9] for more details on the notion of validity of a formula to compute an error bound, additional variants for accurate and efficient error bounds (including one featuring a stabilized EIM), as well as numerical illustrations for a one-dimensional linear diffusion problem and and three-dimensional acoustic scattering problem. The error bound.E<sup>3</sup> can be of particular interest in the following situations: (i) when the stability constant of the original problem is very small (this is the case in many practical problems), (ii) when very accurate solutions are needed, (iii) when considering a nonlinear problem (for which, in some cases, no error bound is possible until a very tight tolerance is reached, see [ 18]).

### **References**


**Open Access** This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

### **Chapter 4 Resources: Software and Tutorials**

### **4.1 Mordicus: Reduced-Order Methods Designed for Industrial Usage**

MORDICUS is the French acronym for "Méthodes d'Ordre RéDuIt Conçues pour des Usages induStriels", translated to "reduced-order methods designed for industrial usage". It is the name of a collaborative project that took place from 2018 to 2023, with the objective of developing a standard for a datamodel and basic computational treatment for reduced-order modeling in the French community.

### *4.1.1 Mordicus Project and Consortium*

"MOR\_DICUS" is a project of type FUI (Fonds Unique Interministériel), in which took part a consortium of ten partners of small and large companies, as well as research institutes and universities. More precisely, the partners are Électricité de France (EDF), Safran Group, ESI Group, Hexagon, Phimeca, Transvalor, CT Ingénierie, CEMOSIS, Armines and Laboratoire Jacques-Louis Lions (LJLL-UPMC), see Fig. 4.1.

The founding is provided by BPI France, pôle systematic and pôle Alsace énergieVie, see Fig. 4.2.

### *4.1.2 Mordicus Library*

Mordicus consists in a Python library and a C++ application, where a datamodel and common basic algorithms and tools are proposed to numerically carry out surrogate modeling and projection-based ROM.

**Fig. 4.2** Funders of the MOR\_DICUS consortium

The sources of the library are available on GitLab.com 1 and the user's license is the GNU LGPLv3. A website containing an installation procedure, a description of the library, some tutorial and a documentation is also available. 2 The rest of the section is inspired from the "Numerical methods" 3 section of this documentation.

The Mordicus library is constructed with the datamodel at its center, and the methods and algorithms act as functions that can modify the state of the datamodel. As detailed in Sect. 2.3, Reduced-Order Model workflows involve an offline (or learning) stage, and an online (or exploitation) stage. The different steps can be vastly adapted and mixed together to produce any complex workflows, as long as the algorithm makes sense. For instance, a regressor can be trained directly on highdimensional data, and the data compression step would not exist. Another example is the Reduced-Basis method, where the notions of offline and online stages merge together, since HFM and ROM resolutions are alternatively computed when constructing the reduced-order basis.

Mordicus provides a datamodel, standard interfaces and basic algorithms, as illustrated in Fig. 4.3. Based on Mordicus, applicative modules can then be developed and propose new workflows and algorithms. The library genericROM, described in Sect. 4.2, is an applicative module of Mordicus.

Simple algorithms are proposed, and can be used by any applicative module.

• SVD: for computing truncated singular value decomposition of lower triangular matrices.

<sup>1</sup> https://gitlab.com/mor\_dicus/mordicus.

<sup>2</sup> https://mordicus.readthedocs.io.

<sup>3</sup> https://mordicus.readthedocs.io/en/latest/\_methods/index.html.

**Fig. 4.3** Organization and overview of the datamodel of the Mordicus library


More details on the datamodel are provided in the next Section.

### *4.1.3 Mordicus Datamodel*

The main feature of Mordicus is a datamodel adapted to reduced-order modeling. It has been constructed to facilitate collaboration, by proposing three main classes CollectionProblemData, ProblemData and Solution, supposed to be populated and handle in the same fashion by all users (see Fig. 4.3):


These classes also contain numerous functions to iterate, modify and handle the data structure. We present a few important ones below:

	- DefineVariabilityAxes(): sets the axes of variability, can be strings for nonparametrized variability, or floats,
	- SnapshotsIterator(): returns an iterator over snapshots of solutions of a given name,

Notice that the functions acting on solution objects automatically iterate over all the problemDatas included in the collectionProblemData.

	- UncompressSolution(): uncompress the reducedCoordinates of a solution of a given name, and update to corresponding solution.snapshots,
	- GetLoadingsOfType(): returns all loadings of a specific type in a list,
	- GetParameterAtTime(): returns the parameter value at a specitiy time (with time interpolation if needed).
	- UncompressSnapshotAtTime(): uncompress the reducedCoordinates of the solution at a given time, and update to corresponding snapshots,
	- GetTimeSequenceFromSnapshots(): returns the time sequence from the snapshots dictionary,
	- ConvertReducedCoordinatesReducedOrderBasisAtTime(): converts the reducedSnapshot at a given time from the current reducedOrder-Basis to a newReducedOrderBasis using a projectedReducedOrderBasis.

The data-model has also been thought to be agile and customizable, by allowing developers to propose other classes, in their applicative module, or variant classes of the ones contained in the subfolders of Containers (e.g. Loadings, Meshes).

### **4.2 GenericROM Library**

The genericROM software consists in a Python library, acting as an applicative module of Mordicus, and developed at Safran.

The sources of the library are available on GitLab.com 4 and the user's license is BSD 3-Clause, a very permissive license that allows the users to use and redistribute the sources, with or without modification. A website containing an installation process, a description of the library, some tutorial and a documentation are also available. <sup>5</sup>

<sup>4</sup> https://gitlab.com/drti/genericrom.

<sup>5</sup> https://genericrom.readthedocs.io.

### *4.2.1 Main Available Methods*

We refer to Sect. 2.3 for an description of the main steps of projection-based reducedorder modeling methods.

The simplest data compression method available in genericROM is the snapshot-POD, as described in Sect. 2.3.3.2. An implementation is available in parallel with distributed memory, by partitioning the domain .Ω in subdomains, as well as an incremental version.

The only hyper-reduction method available to date in genericROM is the ECM, see Sect 2.3.5. A variant of th Nonnegative Orthogonal Matching Pursuit Algorithm 2.2 is implemented, where randomly chosen quadrature point can be added at each iteration to prevent the algorithm to fall into local minima.

The physical problems available for reduction by genericROM are


### *4.2.2 Noninstrusivity and Nonparametrized Variability*

In genericROM, projection-based ROMs can be constructed even when the snapshots are generated by commercial software, as long as readers for the meshes and computed snapshots are available (or can be developed). The handling of projectionbased ROM workflows without having to modify the assembling routines of the reference HFM code is the reason why the library is coined nonintrusive. To do so, the assembling routine of operators and right-hand sides are handling directly in genericROM, using the finite element engine of BasicTools <sup>6</sup> developped at Safran.

This nonintrusive feature is rare in the literature, and has important advantages in an industrial context: expensive and already computed databases of snapshots can be used to construct a ROM. Besides, the ROM library can be updated without having to undergo complex certification of the evolution of reference HFM codes. Moreover, when these codes are commercial, intrusive implementations of projection-based ROM are not possible, or may require expensive contracts with code editors. An example of a ROM constructed with genericROM for a nonlinear transient thermal problem using snapshots computed by Abaqus is available in [ 2].

GenericROM can also deal with nonparametrized variability, in the sense that the variability at the origin of the differences in the snapshots computed in the data

<sup>6</sup> https://gitlab.com/drti/basic-tools.

generation step is not required to be formalized, or even known. This variability takes usually the form of parameters in the definition of the reference HFM, but can also be unknown when for instance the boundary conditions come from a first numerical simulation. Besides, one may want to compute a reduced model by imposing a variability outside a parametrization used to generate the database of snapshots. The nonparametrized variability feature comes with an efficiency trade-off, since, for a certain variability to be changed in a nonparametrized fashion, the ROM needs to assemble the uncompressed version of the corresponding terms in the weak form, before compressing them in the reduced problem (2.13). GenericROM offers the possibility to precompute all the parametrized variability, in order for the ROM to compute efficiently the corresponding terms in the reduced problem, with the consequence that such terms can only be updated following this parametrization. More details on the efficiency are given in the next section. The next two sections are inspired from the "Numerical methods" 7 and "Tutorials" 8 section of the genericROM documentation.

### *4.2.3 Precomputations for Efficiency*

In ROM workflows, the common measure of efficiency is the speedup, defined as the ratio between the computation duration of the HFM and the one of the ROM. For a given accuracy, the higher the speedup, the better. Various elements can be taken into account, including code optimization of the online stage. The methods enabling algorithmic complexity gains have been presented in Sect. 2.3.3. In this section, we give additional implementation elements and illustrations to give some intuition on the mechanisms at play.

Depending on the considered equations and nature of variability, the goal is to precompute as many quantities as possible in the offline stage, to leave as few operations as possible to the online stage. In the HFM, the assembling of the global tangent operator at each Newton iteration, using the high-fidelity quadrature, reads:

$$\frac{D\overline{\mathcal{T}\_{\mu}}}{Du}\left(\boldsymbol{u}^{k}\right)\_{ij} = \sum\_{\mathbf{g}=1}^{N\_{\mathfrak{g}}} \boldsymbol{\omega}\_{\mathfrak{g}}\left[\boldsymbol{\epsilon}\left(\boldsymbol{\varphi}\_{j}\right) : \mathcal{K}\left(\boldsymbol{\epsilon}\left(\boldsymbol{u}^{k}\_{\mu}\right), \mathbf{y}\_{\mu}\right) : \boldsymbol{\epsilon}\left(\boldsymbol{\varphi}\_{i}\right)\right] \left(\mathbf{x}\_{\mathfrak{g}}\right), \ 1 \le i, \ j \le N. \ (4.1)$$

In the ROMs implemented in genericROM, the assembling of the reduced global tangent operator at each Newton iteration, using the reduced quadrature computed by ECM, reads:

.

<sup>7</sup> https://genericrom.readthedocs.io/en/latest/\_methods/index.html.

<sup>8</sup> https://genericrom.readthedocs.io/en/latest/\_tutorials/Mechanical/MecaSequential/index.html.

$$\frac{D\overline{\mathcal{F}\_{\mu}}}{D\hat{u}}\left(\hat{u}^{k}\right)\_{ij} \approx \sum\_{\mathbf{g}=1}^{n\_{\mathbf{g}}} \hat{\boldsymbol{\omega}}\_{\mathbf{g}}\left[\boldsymbol{\epsilon}\left(\boldsymbol{\psi}\_{j}\right) : \mathcal{K}\left(\boldsymbol{\epsilon}(\hat{\boldsymbol{u}}\_{\mu}^{k}), \mathbf{y}\_{\mu}\right) : \boldsymbol{\epsilon}\left(\boldsymbol{\psi}\_{i}\right)\right] \left(\hat{\boldsymbol{x}}\_{\mathbf{g}}\right), \ 1 \le i, j \le n \ll N. \tag{4.2}$$

The corresponding linear systems are illustrated on Figs. 4.4 and 4.5. In particular, . *n* should be much smaller than.*N* to obtain a significant speedup despite the fact that the linear system is sparse in the HFM and dense in the ROM, see Sect. 2.3.6.

Based on Figs. 4.4 and 4.5, the speedup is clear for the resolution part of the linear system. We still need to illustrate how precomputations in the offline phase allow an efficient assembling of the reduced problem. Denote . *d* the number of unknowns of second-order tensor dual quantities. In 3D,.*d* = 6, for instance for the strain tensor, the unknowns are .Ω11, Ω22, Ω33, Ω12, Ω23, Ω31. The reduced global tangent operator at each Newton iteration, using the reduced quadrature, can be written

$$\begin{aligned} \frac{D\overline{\mathcal{F}}\_{\mu}}{D\hat{\boldsymbol{u}}} \left(\hat{\boldsymbol{u}}^{k}\right)\_{ij} & \approx \\ \sum\_{g=1}^{n\_{\mathcal{S}}} \hat{\boldsymbol{\alpha}}\_{\mathcal{S}} \sum\_{l=1}^{d} \left[\boldsymbol{\epsilon}\_{l} \left(\boldsymbol{\psi}\_{j}\right) \left(\hat{\boldsymbol{x}}\_{\mathcal{S}}\right)\right] \sum\_{m=1}^{d} \left[\mathcal{K}\_{l,m} \left(\boldsymbol{\epsilon}\left(\hat{\boldsymbol{u}}^{k}\_{\mu}\right),\boldsymbol{y}\_{\mu}\right) \left(\hat{\boldsymbol{x}}\_{\mathcal{S}}\right)\right] \left[\boldsymbol{\epsilon}\_{m} \left(\boldsymbol{\psi}\_{i}\right) \left(\hat{\boldsymbol{x}}\_{\mathcal{S}}\right)\right], \; 1 \le i, \; j \le n. \end{aligned} \tag{4.3}$$

In genericROM, the code for computing this quantity in the online stage is:

```
reducedTangentMatrix = np.einsum('g,lgj,glm,mgi->ij', 
    reducedIntegrationWeights, 
    reducedEpsilonAtReducedIntegPoints, localTgtMat, 
    reducedEpsilonAtReducedIntegPoints, optimize = True)
```
where

.

.


The object reducedEpsilonAtReducedIntegPoints is a third-order tensor containing the components of the strain tensor applied to the ROB and evaluated at the reduced quadrature points. This quantity is precomputed in the offline stage. The object localTgtMat is a third-order tensor containing the components of the local tangent matrix evaluated at the reduced quadrature points: this quantity is computed online by the constitutive law solver.

### *4.2.4 Tutorials and Datasets*

Various tutorials detailing some physical uses cases and capabilities of genericROM are available online. 9 Namely, variants of the POD-ECM methods are presented for nonlinear structural mechanics and nonlinear transient thermal problems, featuring various boundary conditions. The datasets needed to run these tutorials are also provided.

We detail here one tutorial, 10 corresponding to a nonlinear structural mechanics case, as described in [ 1]. To run this tutorial, a Z-set <sup>11</sup> license is required for the constitutive law solver (Z-mat).

### **4.2.4.1 Description of the Physics Problem**

Consider a 1 m-wide cube, illustrated in Fig. 4.6. This cube is subjected to


The modeling hypothesis are: quasistatic equilibrium equations for the deformable solid in small perturbations. The material is elastoviscoplastic with nonlinear hardening and a Norton flow. The constitutive law is specified by the .mat file.

<sup>9</sup> https://genericrom.readthedocs.io/en/latest/\_tutorials/index.html.

<sup>10</sup> https://genericrom.readthedocs.io/en/latest/\_tutorials/Mechanical/MecaSequential/index.html.

<sup>11</sup> http://www.zset-software.com.

**Fig. 4.6** (left) Dual mesh of the test case with the accumulated plasticity filed at the end of the simulation, (right) mesh and maximal temperature loading field (genericROM documentation)

### **4.2.4.2 Reduction Strategy**

This is a simple "reproducting use case", where a reduced-order model is constructed to reproduce a high-fidelity solution and check its quality. This use case contains the modeling complexity of the high-pressure turbine blade cyclic extrapolation, except for the parallel computation in distributed memory and the fact that the turbine features two material (elastic and elastoviscoplastic). The number of simulation cycles can be seen as the input, while the outputs are the displacement field "U" and the accumulated plastic strain field. *p* (named "evrcum").

The reduced-order basis is generated by the snapshot-POD method "Compress-Data" from the "DataCompressors.FusedSnapshotPOD" module.

The operator compression step is computed using the ECM (Empirical Cubature Method) "CompressOperator" from the "OperatorCompressors.Mechanical" module.

### **4.2.4.3 Algorithm**

.

Snapshots are generated by Z-set and read from the Z-set format. The quality of the data compression is evaluated by computed the relative projection error of the high-fidelity solutions on the reduced-order basis:

$$\frac{\left\|\boldsymbol{u} - \sum\_{i=1}^{N} \left(\boldsymbol{u}, \Psi\_{i}^{\boldsymbol{u}}\right) \Psi\_{i}^{\boldsymbol{u}}\right\|}{\left\|\boldsymbol{u}\right\|} \leq 1 \times 10^{-9},\tag{4.4}$$

where.Ψ*<sup>u</sup> <sup>i</sup>* denote the POD modes (or vectors of the reduced-order basis). The same quantity is considered from. *p*.

The quality of the complete procedure is evaluated by compute the relative .l2 norm error between the high-fidelity solutions and the reduced ones:

$$\frac{\left\|\boldsymbol{u} - \sum\_{i=1}^{N} \boldsymbol{\chi}\_{i}^{\boldsymbol{u}} \boldsymbol{\Psi}\_{i}^{\boldsymbol{u}}\right\|}{\left\|\boldsymbol{u}\right\|} \leq 1 \times 10^{-9},\tag{4.5}$$

with.γ *<sup>u</sup> <sup>i</sup>* the coefficients of the reduced solution (or "reducedCoordinates") computed by the reduced-order model.

### **4.2.4.4 Code—Offline Stage**

List of imports required for the offline stage of this example:

.

```
from genericROM.IO import ZsetMeshReader as ZMR 
from genericROM.IO import ZsetSolutionReader as ZSR 
from Mordicus.Containers import ProblemData as PD 
from Mordicus.Containers import CollectionProblemData as CPD 
from Mordicus.Containers import Solution as S 
from genericROM.FE import FETools as FT 
from genericROM.DataCompressors import FusedSnapshotPOD as SP 
from genericROM.OperatorCompressors import Mechanical 
from Mordicus.IO import StateIO as SIO 
import numpy as np
```
Then, filename and dimensions related to the mesh and the solutions have to be declared, readers are initalized, and the mesh is read:

```
folder = 
    GetTestDataPath()+"Zset"+os.sep+"MecaSequential"+os.sep 
meshFileName = folder + "cube.geof" 
solutionFileName = folder + "cube.ut" 
meshReader = ZMR.ZsetMeshReader(meshFileName) 
solutionReader = ZSR.ZsetSolutionReader(solutionFileName) 
mesh = meshReader.ReadMesh() 
numberOfNodes = mesh.GetNumberOfNodes() 
numberOfIntegrationPoints = 
    FT.ComputeNumberOfIntegrationPoints(mesh) 
nbeOfComponentsPrimal = 3 
nbeOfComponentsDual = 6
```
Then, the part of the ECM algorithm depending only on the mesh is carried out:

```
operatorPreCompressionData = 
    Mechanical.PreCompressOperator(mesh)
```
Then, the objects "Solution" are declared and populated with data from the precomputed Z-set solutions:

```
outputTimeSequence = solutionReader.ReadTimeSequenceFromSolutionFile() 
solutionU = S.Solution("U", nbeOfComponentsPrimal, numberOfNodes, primality 
    = True)
```

```
solutionSigma = S.Solution("sigma", nbeOfComponentsDual, 
    numberOfIntegrationPoints, primality = False) 
solutionEvrcum = S.Solution("evrcum", 1, numberOfIntegrationPoints, 
    primality = False) 
for time in outputTimeSequence: 
  solutionU.AddSnapshot(solutionReader.ReadSnapshot("U", time, 
       nbeOfComponentsPrimal, primality=True), time) 
  solutionSigma.AddSnapshot(solutionReader.ReadSnapshot("sig", time, 
       nbeOfComponentsDual, primality=False), time) 
  solutionEvrcum.AddSnapshot(solutionReader.ReadSnapshotComponent("evrcum", 
       time, primality=False), time)
```
Then, the objects "CollectionProblemData" and "ProblemData" are declared, which will enable to agregate the "Solution" objects previously constructed in a standard fashion in Mordicus:

```
problemData = PD.ProblemData("MecaSequential") 
problemData.AddSolution(solutionU) 
problemData.AddSolution(solutionSigma) 
problemData.AddSolution(solutionEvrcum) 
collectionProblemData = CPD.CollectionProblemData() 
collectionProblemData.AddVariabilityAxis('config', str, 
    description="dummy variability") 
collectionProblemData.DefineQuantity("U", "displacement", "m") 
collectionProblemData.DefineQuantity("sigma", "stress", "Pa") 
collectionProblemData.DefineQuantity("evrcum", "accumulated 
    plasticity", "") 
collectionProblemData.AddProblemData(problemData, 
    config="case-1")
```
Then, the.*L*2(Ω) correlation operator between snapshots is computed (identified by "U"):

```
snapshotCorrelationOperator = 
    {"U":FT.ComputeL2ScalarProducMatrix(mesh, 3)}
```
Then, using the snapshots-POD method, we compute the reduced-order basis for the solutions "U" with the.*L*2(Ω) correlation operator, and for the solutions. *p* without correlation operator (the default operator is the identity) as a precomputing step for the Gappy-POD reconstruction method on. *p*.

```
SP.CompressData(collectionProblemData, "U", 1.e-6, 
    snapshotCorrelationOperator["U"]) 
SP.CompressData(collectionProblemData, "evrcum", 1.e-6)
```
Then, we compute the reduced coefficients (or "reducedCoordinates") by projecting the high-fidelity snapshots onto the reduced-order basis:

```
collectionProblemData.CompressSolutions("U", 
   snapshotCorrelationOperator["U"])
```
Notice that the two previous steps can be done in one by setting the attribute "compressSolutions = True" in the function "SP.CompressData". Then, we verify the quality of the data compression on "U":

```
reducedOrderBasisU = collectionProblemData.GetReducedOrderBasis("U") 
CompressedSolutionU = solutionU.GetReducedCoordinates() 
compressionErrors = [] 
for t in outputTimeSequence: 
  reconstructedCompressedSolution = np.dot(CompressedSolutionU[t], 
       reducedOrderBasisU) 
  exactSolution = solutionU.GetSnapshot(t) 
  norml2ExSol = np.linalg.norm(exactSolution) 
  if norml2ExSol != 0: 
     relError = 
          np.linalg.norm(reconstructedCompressedSolution-exactSolution) 
     /norml2ExSol 
  else: 
     relError = 
          np.linalg.norm(reconstructedCompressedSolution-exactSolution) 
     compressionErrors.append(relError)
```
Then, we carry out the ECM algorithm to determine the reduced quadrature scheme:

```
Mechanical.CompressOperator(collectionProblemData, 
    operatorPreCompressionData, mesh, 1.e-5, 
listNameDualVarOutput = ["evrcum"], 
    listNameDualVarGappyIndicesforECM = ["evrcum"])
```
Finally, at the end of the offline, the Modicus datamodel containing the results of this stage, is saved on disk in order to use it during the online stage.

```
SIO.SaveState("collectionProblemData", collectionProblemData) 
SIO.SaveState("snapshotCorrelationOperator", 
    snapshotCorrelationOperator)
```
### **4.2.4.5 Code—Online Stage**

List of imports required for the offline stage of this example:

```
from genericROM.IO import ZsetInputReader as ZIR 
from genericROM.IO import ZsetMeshReader as ZMR 
from genericROM.IO import ZsetSolutionReader as ZSR 
from genericROM.IO import ZsetSolutionWriter as ZSW 
from Mordicus.Containers import ProblemData as PD 
from Mordicus.Containers import Solution as S 
from genericROM.FE import FETools as FT 
from genericROM.OperatorCompressors import Mechanical as Meca 
from Mordicus.IO import StateIO as SIO 
import numpy as np
```
First, data saved on disk at the end of the offline stage is read:

```
collectionProblemData = SIO.LoadState("collectionProblemData") 
operatorCompressionDataMechanical = 
    collectionProblemData.GetOperatorCompressionData("U") 
snapshotCorrelationOperator = 
    SIO.LoadState("snapshotCorrelationOperator") 
reducedOrderBases = 
    collectionProblemData.GetReducedOrderBases()
```
Then, filename and dimensions related to the mesh and the solutions have to be declared and readers are initialized, in the same fashion as the offline stage. We mention here that the physical setting for the online stage has been taken identical to the one used in the offline stage (the folder "MecaSequential/"). In meaningful workflow, the physical setting for the online stage would be different.

```
folder = 
    GetTestDataPath()+"Zset"+os.sep+"MecaSequential"+os.sep 
inputFileName = folder + "cube.inp" 
inputReader = ZIR.ZsetInputReader(inputFileName) 
meshFileName = folder + "cube.geof"
```
Then, the mesh is read (which is required when the variability is not parametrized):

mesh = ZMR.ReadMesh(meshFileName)

Then, an object "ProblemData" is defined, which will store the data computed during the online stage:

```
onlineProblemData = PD.ProblemData("Online") 
onlineProblemData.SetDataFolder(os.path.relpath(folder, 
    folderHandler.scriptFolder))
```
Then, the temporal sequence and the constitutive law are read from the Z-Set input file. These are "fixed data" for the online resolution:

```
timeSequence = inputReader.ReadInputTimeSequence() 
constitutiveLawsList = 
    inputReader.ConstructConstitutiveLawsList() 
onlineProblemData.AddConstitutiveLaw(constitutiveLawsList)
```
Then, the loadings and initial condition are read from the Z-Set input file and are reduced by projecting them onto the reduced-order basis:

```
loadingList = inputReader.ConstructLoadingsList() 
onlineProblemData.AddLoading(loadingList) 
for loading in onlineProblemData.GetLoadingsForSolution("U"): 
  loading.ReduceLoading(mesh, onlineProblemData, 
      reducedOrderBases, operatorCompressionData)
```

```
initialCondition = inputReader.ConstructInitialCondition() 
onlineProblemData.SetInitialCondition(initialCondition) 
initialCondition.ReduceInitialSnapshot(reducedOrderBases, 
    snapshotCorrelationOperator)
```
Then, the reduced solution is computed in a nonintrusive fashion using a reduced Newton iterative algorithm for solving the reduced nonlinear system of equations at each time-step:

```
onlineCompressedSolution = 
    Meca.ComputeOnline(onlineProblemData, timeSequence, 
    operatorCompressionDataMechanical, 1.e-8)
```
Then, the reduced coefficients (or "reducedCoordinates") for the dual quantified of interest, here. *p*, are computed using the online part of the Gappy-POD:

```
onlineData = onlineProblemData.GetOnlineData("U") 
onlineEvrcumCompressedSolution, errorGappy = 
    Meca.ReconstructDualQuantity("evrcum", 
    operatorCompressionDataMechanical, 
onlineData, timeSequence = 
    list(onlineCompressedSolution.keys()))
```
In order to compare the reduced solutions to the high-fidelity reference ones, "Solution" objects are created and populated with precomputed Z-set solutions:

```
numberOfIntegrationPoints = 
    FT.ComputeNumberOfIntegrationPoints(mesh) 
nbeOfComponentsPrimal = 3 
numberOfNodes = mesh.GetNumberOfNodes() 
solutionFileName = folder + "cube.ut" 
solutionReader = ZSR.ZsetSolutionReader(solutionFileName) 
outputTimeSequence = 
    solutionReader.ReadTimeSequenceFromSolutionFile() 
solutionEvrcumExact = S.Solution("evrcum", 1, 
    numberOfIntegrationPoints, primality = False) 
solutionUExact = S.Solution("U", nbeOfComponentsPrimal, 
    numberOfNodes, primality = True) 
for t in outputTimeSequence: 
  evrcum = solutionReader.ReadSnapshotComponent("evrcum", t, 
      primality=False) 
  solutionEvrcumExact.AddSnapshot(evrcum, t) 
  U = solutionReader.ReadSnapshot("U", t, 
      nbeOfComponentsPrimal, primality=True) 
  solutionUExact.AddSnapshot(U, t)
```
"Solution" objects corresponding to the reduced solutions are constructed, populated with the reduced coefficients (or "reducedCoordinates") computed by the online stage, and reconstructed on the complete mesh:

```
solutionEvrcumApprox = S.Solution("evrcum", 1, numberOfIntegrationPoints, 
    primality = False) 
solutionEvrcumApprox.SetReducedCoordinates(onlineEvrcumCompressedSolution) 
solutionEvrcumApprox.UncompressSnapshots(reducedOrderBases["evrcum"]) 
solutionUApprox = S.Solution("U", nbeOfComponentsPrimal, numberOfNodes, 
    primality = True) 
solutionUApprox.SetReducedCoordinates(onlineCompressedSolution) 
solutionUApprox.UncompressSnapshots(reducedOrderBases["U"])
```
Then, we verify the quality of the reduced solutions "U" and "evrcum":

```
ROMErrorsU = [] 
ROMErrorsEvrcum = [] 
for t in outputTimeSequence: 
  exactSolution = solutionEvrcumExact.GetSnapshotAtTime(t) 
  approxSolution = solutionEvrcumApprox.GetSnapshotAtTime(t) 
  norml2ExactSolution = np.linalg.norm(exactSolution) 
  if norml2ExactSolution > 1.e-10: 
     relError = 
          np.linalg.norm(approxSolution-exactSolution)/norml2ExactSolution 
  else: 
     relError = np.linalg.norm(approxSolution-exactSolution) 
  ROMErrorsEvrcum.append(relError) 
  exactSolution = solutionUExact.GetSnapshotAtTime(t) 
  approxSolution = solutionUApprox.GetSnapshotAtTime(t) 
  norml2ExactSolution = np.linalg.norm(exactSolution) 
  if norml2ExactSolution > 1.e-10: 
     relError = 
          np.linalg.norm(approxSolution-exactSolution)/norml2ExactSolution 
  else: 
     relError = np.linalg.norm(approxSolution-exactSolution) 
  ROMErrorsU.append(relError)
```
Finally, reduced predictions for "U" and "evrcum" are exported in the Z-set format:

```
onlineProblemData.AddSolution(solutionUApprox) 
onlineProblemData.AddSolution(solutionEvrcumApprox) 
ZSW.WriteZsetSolution(mesh, meshFileName, "reduced", 
    collectionProblemData, onlineProblemData, "U")
```
### **4.2.4.6 Results**

A comparison between the reduced and reference high-fidelity solutions is illustrated in Fig. 4.7.

**Fig. 4.7** Comparison between the reduced and reference high-fidelity solutions: (top) on the displacement "U", (bottom) on the accumulated plasticity "p" (genericROM documentation)

### **References**


**Open Access** This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

# **Chapter 5 Industrial Application: Uncertainty Quantification in Lifetime Prediction of Turbine Blades**

This chapter is a synthesis of the previous ones, since many introduced concepts are applied herein. The complete ROM-net workflow, described in Sect. 2.4.2 is applied to the quantification of the uncertainty of dual quantities (such as the accumulated plastic strain and the stress tensor) on an real-life turbine blade, generated by the uncertainty of the temperature loading field. The numerical experiments make use of the codes Mordicus and genericROM, introduced respectively in Sects. 4.1 and 4.2. The content of this chapter is inspired from our publication [ 9].

Computing the fatigue lifetime of one such blade requires simulating its behavior until the stabilization of the mechanical response, which last several weeks using Abaqus [ 23] because of the size of the mesh, the complexity of the constitutive equations, and the number of loading cycles in the transient regime. With such a computation time, uncertainty quantification with the Monte Carlo method is unaffordable. In addition, such simulations are too time-consuming to be integrated in design iterations, which limits them to the final verification steps, while the design process still relies on simplified models. Accelerating these complex simulations is a key challenge while maintaining a satisfying accuracy, as it would provide useful numerical tools to improve design processes and quantify the effect of the uncertainties on the environment of the system.

Simulations are accelerated using a dictionary of reduced order models, with a classifier able to select which local reduced order model to be used for a new temperature loading. A dataset of 200 solutions is computed in a Finite Element approximation space of dimension in the order of the million, for various instances of the temperature field loading, in parallel in 7 days and 9 h on 48 cores. These solutions are computed over 11 time steps in the first cycle, using a scalable Adaptive MultiPreconditioned FETI (AMPFETI) solver [ 3] in Z-set finite-element software [ 17]. The dataset is partitioned into two clusters using a k-medoids algorithm with a ROMoriented dissimilarity measure in 5 min; the corresponding local ROMs, using POD data compression and ECM operator compression, are trained in 2 h and 30 min. An automatic reduced model recommendation procedure, allowing to decide which local ROM to be used for a new temperature loading, is trained in the form of a logistic regression classifier in 16 min. A meta-model is used to reconstruct the dual quantities of interest over the complete mesh from their values on the reduced integration points, in the form of a multi-task Lasso, which takes 1 h to train for 14 dual fields. The uncertainties on dual quantities of interest, such as the accumulated plastic strain and the stress tensor, are quantified by using our trained ROM-net on 1008 Monte Carlo draws of the temperature loading field in 2 h and 48 min, which corresponds to a speedup greater than 600 with respect to our highly optimized domain decomposition AMPFETI solver. Expected values for the Von Mises stress and the accumulated plastic strain have 0.99-confidence intervals width of respectively 1.66% and 2.84% (relative to the corresponding prediction for the expected value). As a verification stage, 20 reference solutions are computed on new temperature loadings, and dual quantities of interest are predicted with relative accuracy in the order of 1% to 2%, while the location of the maximum value is perfectly predicted.

In what follows, we describe the industrial dataset, the hypotheses of the model, and the objective of the present study. The proposed workflow for uncertainty quantification is then applied on this industrial configuration.

### **5.1 Industrial Context**

The industrial test case of interest consists in predicting the mechanical behavior of a high-pressure (HP) turbine blade in an aircraft engine with uncertainties on the thermal loading. The industrial context and the models for the mechanical behavior and the thermal loading are presented, with a particular emphasis on the assumptions that have been made. For confidentiality reasons, mesh sizes and numerical values corresponding to the industrial dataset are not given, and the provided figures and plots do not contain any color map or physical numerical value. The accuracy of the predictions made by our methodology are given in the form of relative errors.

### *5.1.1 Thermomechanical Fatigue of High-Pressure Turbine Blades*

High-pressure turbine blades are critical parts in an aircraft engine. Located downstream of the combustion chamber, they are subjected to extreme thermomechanical loadings resulting from the combination of centrifugal forces, pressure loads, and hot turbulent fluid flows whose temperatures are higher than the material's melting point. The repeating thermomechanical loading over time progressively damages the blades and leads to crack initiation under thermomechanical fatigue. Predicting the fatigue lifetime is crucial not only for safety reasons, but also for ecological issues, since reducing fuel consumption and improving the engine's efficiency requires increasing the temperature of the gases leaving the combustion chamber.

High-pressure turbine blades are made of monocrystalline nickel-based superalloys that have good mechanical properties at high temperatures. To reduce the temperature inside this material, the blades contain cooling channels where flows relatively fresh air coming from the compressor. In addition, the blade's outer surface is protected by a thin thermal barrier coating. In spite of these advanced cooling technologies, the rotor blades undergo centrifugal forces at high temperatures, causing inelastic strains. Under this cyclic thermomechanical loading repeated over the flights, the structure has a viscoplastic behavior and reaches a viscoplastic stabilized response, where the dissipated energy per cycle still has a nonzero value. This is called *plastic shakedown*, and leads to *low-cycle fatigue*. At cruise flight, the persistent centrifugal force applied at high temperature induces progressive (or time-dependent) inelastic deformations: this phenomenon is called *creep*. In addition, the difference between gas pressures on the extrados and the intrados of the blade generates bending effects. Environmental factors may also locally modify the chemical composition of the material, leading to its *oxidation*. As oxidized parts are more brittle, they facilitate crack initiation and growth. *Thermal fatigue* resulting from temperature gradients is another life-limiting factor. Temperature gradients make colder parts of the structure prevent the thermal expansion of hotter parts, creating thermal stresses. Due to their higher temperatures, the hot parts are more viscous and have a lower yield stress, which make them prone to develop inelastic strains in compression. When the temperature cools down after landing, tensile *residual stresses* appear in parts which were compressed at high temperatures and favor crack nucleation. Given the complex temperature field resulting from the internal cooling channels and the turbulent gas flow, thermal fatigue has a strong influence on the turbine blade's lifetime. In particular, during transient regimes such as take-off, an important temperature gradient appears between the leading edge and the trailing edge of the blade, since the latter has a low thermal inertia due to its small thickness and thus warms up faster.

In short, the behavior of a high pressure turbine blade results from a complex interaction between low-cycle fatigue, thermal fatigue, creep, and oxidation. Due to the cost and the complexity of experiments on parts of an aircraft engine, numerical simulations play a major role in the design of high-pressure turbine blades and their fatigue lifetime assessment. All this knowledge have been learned by scientist and engineers during last decades. In the proposed approach to machine learning for model order reduction, all this knowledge is preserved in local ROMs. It is even more than that, the uncertainty propagation comes to complete this valuable traditional knowledge. We do not expect from artificial intelligence to learn everything in our modeling process.

### *5.1.2 Industrial Dataset and Objectives*

Figure 5.1 gives the geometry and the finite-element mesh of a real high-pressure turbine blade. The mesh is made of quadratic tetrahedral elements, and contains a number of nodes in the order of the million. The elasto-viscoplastic mechani-

**Fig. 5.1** High-pressure turbine blade geometry and mesh (micro-perforations are not modeled) [ 9]

cal behavior is described by a crystal plasticity model presented in the appendix of [ 9]. As explained above, Monte Carlo simulations using a commercial software as Abaqus are unaffordable. With the help of domain decomposition methods, the computation time can be reduced by solving equilibrium equations in parallel on different subdomains of the geometry. Using the implementation of the Adaptive MultiPreconditioned FETI solver [ 3] in Z-set finite-element software [ 17], the simulation of one single loading cycle of the HP turbine blade with 48 subdomains takes approximately 53 min.

The objective is to use a ROM-net to quantify uncertainties on the mechanical behavior of the high-pressure turbine blade, given uncertainties on the thermal loading. The reduction of the computation time should enable Monte Carlo simulations for uncertainty quantification. In particular, we are not interested in predicting the state of the structure after a large number of flight-representative loading cycles. Only one cycle is simulated. Cyclic extrapolation of the behavior of a high-pressure turbine blade has been studied in [ 4] and is out of the scope of this section.

### *5.1.3 Modeling Assumptions*

It is assumed that the heat produced or dissipated by mechanical phenomena has negligible effects in comparison with thermal conduction, which enables avoiding

strongly coupled thermomechanical simulations and running thermal and mechanical simulations separately instead. Under a weak thermomechanical coupling, the first step consists in solving the heat equation to determine the temperature field and its evolution over time. The temperature field history defines the thermal loading and is used to compute thermal strains and temperature-dependent material parameters for the mechanical constitutive laws. Once the thermal loading is known, the temperature-dependent mechanical problem must be solved in order to predict the mechanical response of the structure.

The thermomechanical loading applied to the high-pressure turbine blade during its whole life is modeled as a cyclic loading, with one cycle being equivalent to one flight. The rotation speed of the turbine's rotor is proportional to a periodic function of time.ω(*t*) whose evolution over one period (or cycle, see Fig. 5.2) is representative of one flight with its three main regimes, namely take-off, cruise, and landing. The period (or duration of one cycle) is denoted by. *tc*. The rotation speed between flights . *k* and .*k* + 1 is zero, which means that .ω(*ktc*) = 0 for any integer . *k*. The rotation speed.ω(*t*) is scaled so that its maximum is. 1.

Let .<sup>Ω</sup> <sup>⊂</sup> <sup>R</sup><sup>3</sup> denote the solid body representing the high-pressure turbine blade, with .∂Ω denoting its outer surface. Let .∂Ω*<sup>p</sup>* ⊂ ∂Ω be the surface corresponding to the intrados and extrados. The thermal loading is defined as:

$$\forall \xi \in \Omega, \quad \forall t \in \mathbb{R}\_{+}, \qquad T(\xi, t) = (1 - \omega(t))T\_0 + \omega(t)T\_{\text{max}}(\xi), \tag{5.1}$$

where.*T*<sup>0</sup> = 293 K and.*T*max is the temperature field obtained when the rotation speed reaches its maximum. This field.*T*max is obtained either by an aerothermal simulation or by a stochastic model, as explained later. Similarly, the pressure load applied on .∂Ω*<sup>p</sup>* reads:

$$\forall \xi \in \partial \Omega^p, \quad \forall t \in \mathbb{R}\_+, \qquad p^{\partial \Omega}(\xi, t) = (1 - \omega(t)) p\_0^{\partial \Omega} + \omega(t) p\_{\text{max}}^{\partial \Omega}(\xi), \tag{5.2}$$

where .*p*∂Ω <sup>0</sup> = 1 atm is the atmospheric pressure at sea level, and where .*p*∂Ω max is the pressure field obtained when the rotation speed reaches its maximum. The clamping of the blade's fir-tree foot on the rotor disk is modeled by displacements boundary conditions that are not detailed here.

Small geometric details of the structure have been removed to simplify the geometry. Nonetheless, the main cooling channels are considered. The effects of the thermal barrier coating (TBC) have been integrated in aerothermal simulations, but the TBC is not considered in the mechanical simulation although its damage locally increases the temperature in the nickel-based superalloy and thus affects the fatigue resistance of the structure. Additional centrifugal effects due to the TBC are not taken into account.

The predicted mechanical response of the structure depends on many different factors. Below is a nonexhaustive list of influential factors that are possible sources of uncertainties in the numerical simulation:


For simplification purposes, the only source of uncertainty that is considered in this work is the thermal loading. The equations of the mechanical problem are then seen as parametrized equations, where the parameter is the temperature field .*T*max (see Eq. (5.1)) obtained when the rotation speed reaches its maximum value. The dimension of the parameter space is then the number of nodes in the finiteelement mesh. The mechanical loading is assumed to be deterministic. With the crystal orientation, the constitutive laws and their parameters (or coefficients), they are considered as known data describing the context of the study and given by experts. For details on the constitutive law model, we refer to the appendix of [ 9].

### *5.1.4 Stochastic Model for the Thermal Loading*

A stochastic model is required to take into account the uncertainties on the thermal loading. Given the definition of the thermal loading in Eq. (5.1), we only need to model uncertainties in space through the field.*T*max obtained when the rotation speed reaches its maximum value. The random temperature fields must satisfy some constraints: they must satisfy the heat equation, and they must not take values out of the interval.[0 *K*; *T*melt], where.*T*melt is the melting point of the superalloy. These random fields are obtained by adding random fluctuations to a reference temperature field, see Fig. 5.3. The reference field and comes from aerothermal simulations run with the software *Ansys Fluent*. <sup>1</sup> The data-generating distribution is defined as a Gaussian mixture model made of two Gaussian distributions with the same covariance function but with distinct means, and with a prior probability of .0.5 for each Gaussian distribution. The Gaussian distributions are obtained by taking the four first eigenfunctions of the covariance function (see Karhunen-Loève expansion [ 16]), with a standard deviation of.15 K. Therefore, realizations of the random temperature field read:

$$\forall \xi \in \Omega, \quad T(\xi) = T\_{\text{ref}}(\xi) + \Upsilon\_0 \,\delta T\_0(\xi) + \sum\_{i=1}^{4} \Upsilon\_i \,\delta T\_i(\xi), \tag{5.3}$$

where .*T*ref is the reference field, .δ*T*<sup>0</sup> is a temperature perturbation at the trailing edge whose maximum value is.50 K,.{δ*Ti*}1≤*i*≤<sup>4</sup> are fluctuation modes,.ϒ<sup>0</sup> is a random variable following the Bernoulli distribution with parameter.0.5, and. {ϒ*i*}1≤*i*≤<sup>4</sup> are independent and identically distributed random variables following the standard normal distribution .N(0, 1). The variable .ϒ<sup>0</sup> is also independent of the other variables .ϒ*<sup>i</sup>* . The different fields involved in Eq. (5.3) can be visualized in Fig. 5.3. Equation (5.3) defines a mixture distribution with two Gaussian distributions whose means are .*T*ref and .*T*ref + δ*T*0. We voluntarily define this mixture distribution with .δ*T*<sup>0</sup> adding.50 K in a critical zone of the turbine blade in order to check that our cluster analysis can successfully detect two relevant clusters, i.e., one for fields obtained with.ϒ0(θ ) = 0 and one for fields obtained with.ϒ0(θ ) = 1. Indeed, the temperature perturbation.δ*T*<sup>0</sup> is expected to significantly modify the mechanical response of the high-pressure turbine blade. All the fields .{δ*Ti*}<sup>0</sup>≤*i*≤<sup>4</sup> satisfy the steady heat equation like .*T*ref, which ensures that the random fields always satisfy the heat equation under the assumption of a linear thermal behavior. For nonlinear thermal behaviors, Eq. (5.3) would define surface temperature fields that would be used as Dirichlet boundary conditions for the computation of bulk temperature fields. The assumption of a linear thermal behavior is adopted here to avoid solving the heat equation for every realization of the random temperature field.

Let us now give more details about the construction of the fluctuation modes .{δ*Ti*}<sup>1</sup>≤*i*≤4. First, surface fluctuation modes are computed on the boundary.∂Ω using the method given in [ 21] for the construction of random fields on a curved surface.

<sup>1</sup> https://www.ansys.com/products/fluids/ansys-fluent.

**Fig. 5.3** Reference temperature field (on the left), temperature perturbation at the trailing edge (field 0 =.δ*T*0), and fluctuation modes (fields 1 to 4). The fluctuations in the fourth mode are located inside the blade, in the cooling channels [ 9]

The correlation function is defined as a function of the geodesic distance .*dG* along the surface.∂Ω:

$$\rho(\xi, \xi') = \exp\left(-\frac{d\_G(\xi, \xi')}{d\_G^0}\right),\tag{5.4}$$

where.*d*<sup>0</sup> *<sup>G</sup>* is a correlation length. Geodesic distances are computed using the algorithm described in [ 18, 22] and implemented in the Python library *gdist*. <sup>2</sup> A covariance matrix is built by evaluating the correlation function on pairs of nodes of the outer surface of the finite-element mesh, and multiplying the correlation by the constant variance. The four surface modes are then obtained by finding the four eigenvectors corresponding to the largest eigenvalues of the covariance matrix. The steady heat equation with Dirichlet boundary conditions is solved for each of these surface modes to derive the 3D fluctuation modes, using *Z-set* [ 17] finite-element solver. The Python library *BasicTools* <sup>3</sup> developed by SafranTech is used to read the finite-element mesh and write the temperature fields in a format that can be used for simulations on *Z-set*.

<sup>2</sup> https://pypi.org/project/gdist/.

<sup>3</sup> https://gitlab.com/drti/basic-tools.

### **5.2 ROM-net Based Uncertainty Quantification Applied to an Industrial High-Pressure Turbine Blade**

This section develops the different stages of the ROM-net for the industrial test case presented in the previous section. Given our budget of 200 high-fidelity simulations, a dictionary containing two local ROMs is constructed using our clustering procedure. A logistic regression classifier is trained for automatic model recommendation using information identified by feature selection, followed by an alternative to the Gappy-POD for full-field reconstruction of dual quantities. Then, the results of the uncertainty quantification procedure are presented. Finally the accuracy of the ROM-net is validated using simulations for new temperature loadings.

### *5.2.1 Design of Numerical Experiments*

Given the computational cost of high-fidelity mechanical simulations of the highpressure turbine blade, the training data are sampled from the stochastic model for the thermal loading using a design of experiments. Our computational budget corresponds to 200 high-fidelity simulations, so a database of 200 temperature fields must be built. This database includes two separate datasets coming from two independent DoEs:


These DoEs are built with the platform *Lagun*. <sup>4</sup> The fact that these two datasets come from two separate DoEs is beneficial: as each of them is supposed to have good space-filling properties, they are both representative of the possible thermal loading and can therefore be used to define a training set and a test set for a given learning task. For instance, the classifier trained on the Sobol' DoE can be tested

<sup>4</sup> https://gitlab.com/drti/lagun.

**Fig. 5.4** Visualization of the MaxProj LHS DoE. The marginal distributions are represented on the diagonal. The 5D DoE is projected on 2D subspaces for visualization purposes, in order to check space-filling properties in 2D [ 9]

on the MaxProj LHS DoE. The local ROMs built from snapshots belonging to the MaxProj LHS DoE can make predictions on the Sobol' DoE that will be used for the training of the Gappy surrogates, which is relevant since the Gappy surrogates are supposed to analyze ROM predictions on new unseen data in the exploitation phase.

Drawing random temperature fields as defined in Eq. (5.3) requires sampling data from the random variables .{ϒ*i*}<sup>0</sup>≤*i*≤4, where .ϒ<sup>0</sup> follows the Bernoulli distribution with parameter .0.5 and the variables .ϒ*<sup>i</sup>* for .*i* ∈ [[1; 4]] are independent standard normal variables and independent of.ϒ0. Both DoE methods (Maximum Projection LHS and Sobol' sequence) generate point clouds with a uniform distribution in the unit hypercube. Figures 5.4 and 5.5 show the projections onto 2-dimensional subspaces of the 5D point clouds used to build our datasets. The marginal distributions are plotted to check that they well approximate the uniform distribution. These point

**Fig. 5.5** Visualization of the Sobol' DoE. The marginal distributions are represented on the diagonal. The 5D DoE is projected on 2D subspaces for visualization purposes, in order to check space-filling properties in 2D [ 9]

clouds, considered as samples of a random vector.(χ0, χ1, χ2, χ3, χ4) following the uniform distribution on the unit hypercube, are transformed into realizations of the random vector.(ϒ0, ϒ1, ϒ2, ϒ3, ϒ4) using the following transformations:

$$
\Upsilon\_0 = \mathbb{1}\_{\chi\_0 \ge 1/2} \quad \text{and} \quad \forall i \in \llbracket 1; 4 \}, \quad \Upsilon\_i = F^{-1}(\chi\_i), \tag{5.5}
$$

where .*F*−<sup>1</sup> is the inverse of the cumulative distribution function of the standard normal distribution. The resulting samples define the MaxProj dataset and the Sobol' dataset of random temperature fields, using Eq. (5.3). Each temperature field defines a thermal loading, using Eq. (5.1). The 200 corresponding mechanical problems are solved for one loading cycle with the finite-element software *Z-set* [ 17] with the domain decomposition method described in [ 3], with 48 subdomains. The average computation time for one simulation is 53 min.

### *5.2.2 ROM Dictionary Construction*

The 80 simulations associated to the MaxProj dataset are used as clustering data. Loading all the simulation data and computing the pairwise ROM-oriented dissimilarities takes about 5 min. The ROM-oriented dissimilarity defined in [ 7, Definition 3.11] is computed with.*n* = 1, i.e., each simulation is represented by one field. The dataset is partitioned into two clusters using our implementation of PAM [ 14, 15] k-medoids algorithm, with 10 different random initializations for the medoids. The clustering results can be visualized using Multidimensional Scaling (MDS) [ 2]. MDS is an information visualization method which consists in finding a low-dimensional dataset.*Z*<sup>0</sup> whose matrix of Euclidean distances.*d*(*Z*0)is an approximation of the true dissimilarity matrix. *δ*. To that end, a cost function called stress function is minimized with respect to. *Z*:

$$\mathbf{Z}\_0 = \underset{\mathbf{Z}}{\text{arg min}} \ (\boldsymbol{\zeta}(\mathbf{Z}; \boldsymbol{\delta})) = \underset{\mathbf{Z}}{\text{arg min}} \left( \sum\_{i$$

This minimization problem is solved with the algorithm Scaling by MAjorizing a COmplicated Function (SMACOF, [ 10]) implemented in Scikit-Learn [ 19]. Figure 5.6 show the clusters on the MDS representations with the goal-oriented variants of the ROM-oriented dissimilarity measure, applied on the accumulated plastic strain .*p<sup>o</sup>* cum. The clustering results is compared with the expected clusters corresponding to.ϒ<sup>0</sup> = 0 and.ϒ<sup>0</sup> = 1, the latter corresponds to the perturbation.δ*T*<sup>0</sup> being activated. The obtained clusters almost correspond to the expected ones, with only 4 points with wrong labels out of 80, which quantifies the ability of the ROM-oriented dissimilarity measure on the accumulated plasticity to infer the correct value of.ϒ0. The medoids of the two clusters are given in Fig. 5.7. Cluster. 0 contains temperature fields for which.ϒ<sup>0</sup> = 1, while cluster. 1 contains fields for which.ϒ<sup>0</sup> = 0. It can be observed that the quantity of interest clearly differs from one cluster to the other, while the differences are hardly visible on the displacement field. The displacement field combines deformations associated to different phenomena (thermal expansion, elastic strains, viscoplastic strains) that are not necessarily related to damage in the structure, which could explain why the quantity of interest .*p<sup>o</sup>* cum seems to be more appropriate for clustering in this example.

The simulations used for the clustering procedure can directly provide snapshots for the construction of the local ROMs. To control the duration of their training, only 20 simulations are selected to provide snapshots for the each local ROMs, which represents half of the clusters' populations. These simulations are selected in a maximin greedy approach starting from the medoid (see [ 8, Algorithm 2, Stage 2]

**Fig. 5.6** MDS representation of the clustering results using the ROM-oriented dissimilarity measure on the quantity of interest .*p<sup>o</sup>* cum (goal-oriented variant). On the left, the colors correspond to the expected clusters. On the right, the colors correspond to the clusters identified by the clustering algorithm. The positions of the labels. 0 and. 1 coincide with the positions of the clusters' medoids. The MDS relative error.ς (*Z*0; *δ*)/ς (**0**; *δ*) is.12% [ 9]

**Fig. 5.7** The 3 fields on the left correspond to the medoid of cluster . 0, and those on the right correspond to the medoid of cluster . 1. The fields in the first and the third columns show the differences between the medoids' temperature fields and the reference temperature field .*T*ref (the scale is truncated for the first field). The second and the fourth columns show the displacement magnitude field. <sup>√</sup>*u*.*<sup>u</sup>* (top) and the quantity of interest.*p<sup>o</sup>* cum (bottom) [ 9]

for a example of maximin selection). Figure 5.8 shows which simulations have been selected for the construction of the local ROMs.

The local ROMs are built following the methodology described in Sect. 2.3. The snapshot-POD and the ECM are done in parallel with shared memory on 24 cores. The tolerance for the snapshot-POD is set to.10−<sup>8</sup> for the displacement field, and to. 10−<sup>4</sup> for dual variables (the quantity of interest.*p<sup>o</sup>* cum and the six components of the stress tensor). The POD bases for the dual variables will be used for their reconstruction with the Gappy surrogates. The tolerance for the ECM is set to .5 × 10−4. Locals ROMs each contain 18 displacement modes and between 8 and 13 modes for stress components. The first and second local ROMs contain respectively 10 and 12 modes for the quantity of interest.*p<sup>o</sup>* cum. The ECM selects 506 (resp. 510) integration points for the reduced-integration domain of ROM . 0 (resp. . 1). Building one local ROM takes approximately 2 h and 30 min.

### *5.2.3 Automatic Model Recommendation*

In this section, a classifier is trained for the automatic model recommendation task. The 120 temperature fields coming from the Sobol' dataset are used as training data for the classifier. Their labels are determined by finding their closest medoid in terms of the ROM-oriented dissimilarity measure. Hence, for each temperature field of the Sobol' dataset, two dissimilarities are computed: one with the medoid of the first cluster, and one with the medoid of the second cluster. Once trained, the classifier can be evaluated on the 80 labelled temperature fields of the MaxProj dataset.

Each temperature field is discretized on the finite-element mesh, which contains in the order of the million nodes. To reduce the dimension of the input space and facil-

itate the training phase of the classifier, we apply the geostatistical mRMR feature selection algorithm described in [ 8, Algorithm 1] on data from the Sobol' dataset. First, 800 pairs of nodes are selected in the mesh, which takes 18 s. The 800 corresponding redundancy terms are computed with Scikit-Learn [ 19] in less than 3 s. Figure 5.9 plots the values of these redundancy terms versus the Euclidean distance between the nodes. We observe that the correlation between the redundancy mutual information terms and the distance between the nodes is poor, with a lot of noise. This can be due to the fact that the random temperature fields have been built using Gaussian random fields on the outer surface with an isotropic correlation function depending on the geodesic distance along the surface rather than the Euclidean distance. Since the turbine blade is a relatively thin structure, two nodes, one on the intrados and another one on the extrados, can be close to each other in the Euclidean distance, but with totally uncorrelated temperature fluctuations because of the large geodesic distance separating them. On the contrary, two points on the same side of the turbine blade can have correlated temperature variations while being separated by a Euclidean distance in the order of the blade's thickness. The length of the mutual information's high-variance regime seems to correspond to the blade's chord, which supports this explanation. The thinness of the turbine blade induces anisotropy in the correlation function of the bulk Gaussian random field defining the thermal loading, which implies an anisotropic behavior of the mutual information according to [ 8, Property 1]. The use of a local temperature perturbation .δ*T*<sup>0</sup> in conjunction with fluctuation modes having larger length scales may also partially explain the large variance of redundancy terms. Nonetheless, it remains clear that redundancy terms are smaller as for large distances. This trend is captured by a kriging metamodel (Gaussian process regression) trained with Scikit-Learn in a few seconds, with a sum-kernel involving the Matérn kernel with parameter.5/2 (to get a continuous and twice differentiable metamodel) and length scale . 1, and a white kernel to estimate the noise level of the signal. The curve of the metamodel is given in Fig. 5.9. Then, for each node of the finite-element mesh, the mutual information with the label vari-


**Table 5.1** Classification results

able is computed. The computations of these relevance terms (in the order of the million terms) are distributed between 280 cores, which gives a total computation time of 15 min. Among these features,.5, 986 features are preselected by discarding those with a relevance mutual information lower than.0.05. The geostatistical mRMR selects.11 features in 42 s. The corresponding nodes in the finite-element mesh can be visualized in Fig. 5.9.

**Remark 5.1** The metamodel for redundancy terms could be improved by defining it as a function of the precomputed geodesic distances along the outer surface rather than the Euclidean distances. Each finite-element node would be associated to its nearest neighbor on the outer surface before computing the approximate mutual information from geodesic distances.

The classifier is trained on the Sobol' dataset, using the values of the temperature fields at the .11 nodes identified by the feature selection algorithm. The classifier is a logistic regression [ 1, 5, 6] with elastic net regularization [ 24] implemented in Scikit-Learn. The two hyperparameters involved in the elastic net regularization are calibrated using 5-fold cross-validation, giving a value of.0.001 for the inverse of the regularization strength, and.0.4 for the weight of the.*L*<sup>1</sup> penalty term (and thus.0.6 for the.*L*<sup>2</sup> penalty term). Due to the.*L*<sup>1</sup> penalty term, the classifier only uses. 5 features among the .11 input features. The classifier's accuracy, evaluated on the MaxProj dataset to use new unseen data, reaches.98.75%. The confusion matrix indicates that .100% of the test examples belonging to class. 0 have been correctly labeled, and that .2.38% of the test examples belonging to class . 1 have been misclassified. Table 5.1 summarizes the values of precision, recall and F1-score on test data.

### *5.2.4 Surrogate Model for Gappy Reconstruction*

When using hyper-reduction, the ROM calls the constitutive equations solver only at the integration points belonging to the reduced-integration domain. It is recalled that the ECM selected 506 (resp. 510) integration points for the reduced-integration domain of ROM . 0 (resp. . 1), and that the finite-element mesh initially contains a number of integration points in the order of the million. Therefore, after a reduced simulation, dual variables defined at integration points are known only at integration points of the reduced-integration domain. To retrieve the full field, the Gappy-POD [ 11] finds the coefficients in the POD basis that minimize the squared error between the reconstructed field and the ROM predictions on the reduced-integration domain. This minimization problem defines the POD coefficients as a linear function of the predicted values on the reduced-integration domain. Although these coefficients are optimal in the least squares sense, they can be biased by the errors made by the ROM. To alleviate this problem, we propose to replace the common Gappy-POD procedure by a metamodel or *Gappy surrogate*. The inputs and the outputs of the Gappy surrogate are the same as for the Gappy-POD: the input is a vector containing the values of a dual variable on the reduced-integration domain, and the output is a vector containing the optimal coefficients in the POD basis. One Gappy surrogate must be built for each dual variable of interest: in our case, 7 surrogate models per cluster are required, namely one for the quantity of interest .*p<sup>o</sup>* cum and one for every component of the Cauchy stress tensor.

The training data for these Gappy surrogates are obtained by running reduced simulations with the local ROMs, using the thermal loadings of the Sobol' dataset. Indeed, the two local ROMs have been built on the MaxProj dataset, therefore thermal loadings of the Sobol' dataset can play the role of test data for the ROMs. For each thermal loading in the Sobol' dataset, the true high-fidelity solution is already known since it has been computed to provide training data for the classifier. In addition, the exact labels for these thermal loadings are known, which means that we know which local ROM to choose for each thermal loading of the Sobol' dataset. Given ROM predictions on the reduced-integration domain, the optimal coefficients in the POD basis are given by the projections of the true prediction made by the highfidelity model (the finite-element model) onto the POD modes. This provides the true outputs for the Sobol' dataset, which can then be used as a training set for the Gappy surrogates.

Given the high-dimensionality of the input data (there are more than 500 integration points in the reduced-integration domains) with respect to the number of training examples (120 examples), a multi-task Lasso metamodel is used. The hyperparameter controlling the regularization strength is optimized by 5-fold cross-validation. Training the 14 Gappy surrogates (7 for each cluster) takes 1 h. The Gappy surrogates select between .8% and .18% of the integration points in the reduced-integration domains, due to the.*L*<sup>1</sup> regularization. The mean cross-validated coefficients of determination are.0.9637 (resp..0.8935) for the quantity of interest for cluster. 0 (resp. cluster. 1), and range from .0.9404 to .0.9938 for stress components. These satisfying results mean that it is not required to train a kriging metamodel with the variables selected by Lasso to get nonlinear Gappy surrogates. The Gappy surrogates are then linear, just as the Gappy-POD.

The accuracy gains provided by the Gappy surrogates with respect to classical Gappy-POD on the present industrial case is investigated in Fig. 5.10. Here, . 24 high-fidelity simulation in the first cluster are computed (with.ϒ<sup>0</sup> = 0) as reference, and Gappy surrogates (using two meta models Lasso and ElasticNet) and classical Gappy-POD are computed using ROM . 0. For both variants of meta models and

**Fig. 5.10** Mean over the complete mesh of dual quantities of interest: accumulated plastic strain .*p<sup>o</sup>* cum (left) and component.33 of the stress tensor.σ<sup>33</sup> (right), plotted as points where the x-coordinate is the reference value, and the y-coordinate is the considered reduced prediction [ 9]

both quantities of interest: accumulated plastic strain and the component .33 of the stress tensor, Gappy surrogates provides more accurate predictions than the classical Gappy-POD.

**Remark 5.2** In this strategy, the local ROMs solve the equations of the mechanical problem, which enables using linear surrogate models to reconstruct dual variables. Using surrogate models from scratch instead of local ROMs would have been more difficult, given the nonlinearities of this mechanical problem and the lack of training data for regression. In addition, such surrogate models would require a parametrization of the input temperature fields, whereas the local ROMs use the exact values of the temperature fields on the RID without assuming any model for the thermal loading.

The dictionary-based ROM-net used for mechanical simulations of the highpressure turbine blade is made of a dictionary of two local hyper-reduced order models and a logistic regression classifier. The classifier analyzes the values of the input temperature field at 11 nodes only, identified by our feature selection strategy. For a given thermal loading in the exploitation phase, after the reduced simulation with the local ROM recommended by the classifier, linear cluster-specific Gappy surrogates reconstruct the full dual fields (quantity of interest and stress components) from their predicted values on the reduced-integration domain.

### *5.2.5 Uncertainty Quantification Results*

Once trained, the ROM-net can be applied for the quantification of uncertainties on the mechanical behavior of the HP turbine blade resulting from the uncertainties on the thermal loading. Since the ROM-net online operations can be performed sequentially on one single core, 24 cores are used in order to compute the solution for 24 thermal loadings at once. This way, 42 batches of 24 Monte Carlo simulations

**Fig. 5.11** Histograms and probability density functions of the quantities of interest.*p<sup>o</sup>* cum (left) and .σeq (right) [ 9]

**Table 5.2** Widths of the confidence intervals (CI) for the expectations, expressed as percentages of the estimated expectations


are run in 2 h and 48 min. The 1008 thermal loadings used for this study are generated by randomly sampling points from the uniform distribution on the 5D unit hypercube and applying the transformation given in Eq. (5.5).

The expected values of .*p<sup>o</sup>* cum and .σeq are estimated with the empirical means .*Zn* <sup>=</sup> <sup>1</sup> *n* ∑*<sup>n</sup> <sup>i</sup>*=<sup>1</sup> *Zi* , where.*Zi* are the corresponding samples. The variances of.*p<sup>o</sup>* cum and .σeq are computed using the unbiased sample variance .*S*<sup>2</sup> *<sup>n</sup>* <sup>=</sup> <sup>1</sup> *n* − 1 ∑*n i*=1 ( *Zi* − *Zn* )2 . The Central Limit Theorem gives asymptotic confidence intervals for the expected values: for all.α ∈]0; 1[,

$$I\_n = \left[ \overline{Z}\_n - \phi\_{1-\frac{\sigma}{2}} \sqrt{S\_n^2/n}; \,\overline{Z}\_n + \phi\_{1-\frac{\sigma}{2}} \sqrt{S\_n^2/n}, \right],\tag{5.7}$$

where.φ*<sup>r</sup>* denotes the quantile of order. *r* of the standard normal distribution.N(0, 1), and .*In* is an asymptotic confidence interval with confidence level .1 − α for the expectation. η:.lim*<sup>n</sup>*→+∞ <sup>P</sup>(η <sup>∈</sup> *In*) <sup>=</sup> <sup>1</sup> <sup>−</sup> <sup>α</sup>. The widths of the confidence intervals are expressed as a percentage of the estimated value for the expectations in Table 5.2.

The probability density functions of the quantities of interest can be estimated using Gaussian kernel density estimation (see Sect. 6.6.1. of [ 12]). Figure 5.11 gives

**Fig. 5.12** Workflow for the ROM-net methology applied to the considered industrial setting [ 9]

the histograms and estimated distributions for.*p<sup>o</sup>* cum and.σeq. The shapes of these distributions highly depend on the assumptions made for the stochastic thermal loading.

### *5.2.6 Workflow*

Figure 5.12 provides an illustration of the workflow and the computational time of each step presented above.

### *5.2.7 Verification*

For verification purposes, the accuracy of the ROM-net is evaluated on 20 Monte Carlo simulations with 20 new thermal loadings. These thermal loadings are generated by randomly sampling points from the uniform distribution on the 5D unit hypercube, and applying the transformation given in Eq. (5.5). The reduced simula-

tions are run on single cores. The total computation time for generating a new thermal loading on the fly, selecting the corresponding reduced model, running one reduced simulation and reconstructing the quantities of interest is 4 min on average. As a comparison, one single high-fidelity simulation with *Z-set* [ 17] with 48 subdomains takes 53 min, which implies that the ROM-net computes.13.25 times faster. However, one high-fidelity simulation requires 48 cores for domain decomposition, whereas the ROM-net works on one single core. Hence, using 48 cores to run 48 reduced simulations in parallel,.636 reduced simulations can be computed in 53 min with the ROM-net, while the high-fidelity model only runs one simulation. In addition to the acceleration of numerical simulations, energy consumption is reduced by a factor of .636 in the exploitation phase. In spite of the fast development of high-performance computing, numerical methods computing approximate solutions at reduced computational resources and time are particularly important for many-query problems such as uncertainty quantification, where the intensive use of computational resources is a major concern. Model order reduction and ROM-nets play a prominent role toward *green* numerical simulations [ 20]. Of course, the number of simulations in the exploitation phase must be large enough to compensate the efforts made in the training phase, like in any machine learning or model order reduction problem.

Figures 5.13 and 5.14 show the results for two simulations belonging to cluster . 0 and cluster . 1 respectively. These figures give the difference between the current temperature field as the reference one, i.e., the field.*T* − *T*ref, and the resulting variations of the quantity of interest predicted by the ROM-net and the high-fidelity model, i.e.,.*po*,ROM cum (*T* ) − *po*,HF cum (*T*ref) and.*po*,HF cum (*T* ) − *po*,HF cum (*T*ref). The signs and the positions of the variations of the quantity of interest seem to be quite well predicted by the ROM-net.

Let us introduce a zone of interest .Ω' defined by all of the integration points at which .*p<sup>o</sup>* cum is higher than .0.4 × max *p<sup>o</sup>* cum(*ξ* ) for the thermal loading defined by .*T*ref + δ*T*0. This zone of interest contains 209 integration points. The values of the variables.*p<sup>o</sup>* cum and.σeq averaged over.Ω' are denoted by.*p<sup>o</sup>* cum and.σeq. Table 5.3 gives different indicators quantifying the errors made by the ROM-net: the .*L*<sup>2</sup> relative errors on the whole domain.Ω and on the zone of interest.Ω' , the.*L*<sup>∞</sup> relative errors on .Ω and .Ω' , the relative errors on .*p<sup>o</sup>* cum and .σeq, and the errors on the locations of the points where the fields .*p<sup>o</sup>* cum and .σeq reach their maxima. All the relative errors remain in the order of.1% or.2%, which validates the methodology. In addition, the ROM-net perfectly predicts the position of the critical points at which.*p<sup>o</sup>* cum and. σeq reach their maxima. Figure 5.15 shows errors on the quantities of interest.

**Fig. 5.13** Comparison between high-fidelity predictions (middle column) and ROM-net's predictions (right-hand column). The field on the left represents the difference between the current temperature field (belonging to cluster . 0) and the reference one. The other fields correspond to the increments of the quantity of interest.*p<sup>o</sup>* cum with respect to its reference state obtained with the reference temperature field [ 9]


**Table 5.3** Error indicators for the evaluation of the ROM-net on 20 new thermal loadings

**Fig. 5.14** Comparison between high-fidelity predictions (middle column) and ROM-net's predictions (right-hand column). The field on the left represents the difference between the current temperature field (belonging to cluster . 1) and the reference one. The other fields correspond to the increments of the quantity of interest.*p<sup>o</sup>* cum with respect to its reference state obtained with the reference temperature field [ 9]

**Fig. 5.15** Errors on the quantity of interest.*p<sup>o</sup>* cum. The red (resp. blue) color is used for zones where the quantity of interest is overestimated (resp. underestimated) [ 9]

### **References**


**Open Access** This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

## **Chapter 6 Applications and Extensions: A Survey of Literature**

This chapter contains a literature survey of the work published by the authors in the timeframe of their collaboration, where the concept presented in this book have been applied to real-life industrial settings, and new methodologies have been developed.

The listed contributions are grouped in the following themes: linear manifold learning in Sect. 6.1, nonlinear dimensionality reduction via auto-encoder in Sect. 6.2, piecewise linear dimensionality reduction via dictionary-based ROM-nets in Sect. 6.3 and manifold learning of physics problems assisted by black-box regressors in Sect. 6.4.

### **6.1 Linear Manifold Learning**

**A priori hyper-reduction method: an adaptive approach** Model reduction methods are usually based on offline preliminar simulations to build the shape functions of the reduced order model (ROM) before the computation of the reduced state variables. They are a posteriori approaches. Most of the time these offline computations are as complex as the simulation which we want to simplify by the ROM. The a priori reduction method proposed in [ 17] avoids such preliminary computations. It is an a priori approach based on the analysis of some state evolutions, such that all the state evolutions needed to perform the model reduction are described by an approximate ROM. The ROM and the state evolution are simultaneously improved by the method, thanks to an adaptive strategy. An initial set of known shape functions can be used to define the ROM to adapt, but it is not necessary. The adaptive procedure includes extensions of the subspace spanned by the shape functions of the ROM and selections of the most relevant shape functions in order to represent the state evolution. The hyper-reduction is achieved by selecting a part of the integration points of the finite element model to forecast the evolution of the reduced state variables. In this method, both the number of degrees of freedom and the number of integration points are reduced.

**A nonintrusive distributed reduced-order modeling framework for nonlinear structural mechanics—Application to elastoviscoplastic computations** In [ 8] is proposed a framework that constructs reduced-order models for nonlinear structural mechanics in a nonintrusive fashion and can handle large-scale simulations. Three steps are carried out: (i) the production of high-fidelity solutions by commercial software, (ii) the offline stage of the model reduction, and (iii) the online stage where the reduced-order model is exploited. The nonintrusivity assumes that only the displacement field solution is known, and the proposed framework carries out operations on these simulation data during the offline phase. The compatibility with a new commercial code only needs the implementation of a routine converting the discretized solution into the used data format. The nonintrusive capabilities of the framework are demonstrated on numerical experiments using commercial versions of Z-set and Ansys Mechanical. The nonlinear constitutive equations are evaluated by using an external plugin. The large-scale simulations are handled using domain decomposition and parallel computing with distributed memory. The features and performances of the framework are evaluated on two numerical applications involving elastoviscoplastic materials: the second one involves a model of high-pressure blade, where the framework is used to extrapolate cyclic loadings in 6.5 h, whereas the reference high-fidelity computation would take 9.5 days.

**Fast computation of transient thermal profiling of high-pressure compressors under non-parametrized boundary conditions variability** In [ 9], a transient thermal problem is considered, with a nonlinear term coming from the radiation boundary condition and a nonparametrized variability in the form complex scenarios for the initial condition and the convection coefficients and external temperatures. A posteriori reduced order modeling by snapshot Proper Orthogonal Decomposition is used. To treat the nonlinearity, hyperreduction is required since precomputing the polynomial nonlinearities becomes too expensive for the radiation term. The Empirical Cubature Method, originally proposed for nonlinear structural mechanics, is derived for these equations. The method is applied to the design of high-pressure compressors for civilian aircraft engines, where a fast evaluation of the solution temperature is required when testing new configurations. When using in the reduced solver the same model as the one from the high-fidelity code, the approximation is very accurate. However, when using a commercial code to generate the high-fidelity data, where the implementation of the model and solver is unknown, the reduced model is less accurate but still within engineering tolerances. Hence, the regularizing property of reduced order models, together with a nonintrusive approach, enables the use of commercial software to generate the data, even under some degree of uncertainty in the proprietary model or solver of the commercial software.

**Time Stable Reduced Order Modeling by an Enhanced Reduced Order Basis of the Turbulent and Incompressible 3D Navier–Stokes Equations** In [ 3], the problem of constructing a time stable reduced order model of the 3D turbulent and incompressible Navier–Stokes equations is considered. The lack of stability associated with the order reduction methods of the Navier–Stokes equations is a wellknown problem and, in general, it is very difficult to account for different scales of a turbulent flow in the same reduced space. To remedy this problem, a new stabilization technique is proposed, based on an a priori enrichment of the classical proper orthogonal decomposition (POD) modes with dissipative modes associated with the gradient of the velocity fields. The main idea is to be able to do an a priori analysis of different modes in order to arrange a POD basis in a different way, which is defined by the enforcement of the energetic dissipative modes within the first orders of the reduced order basis. This enables the modeling of the production and dissipation of the turbulent kinetic energy (TKE) in a separate fashion within the high ranked new velocity modes, hence to ensure good stability of the reduced order model. The importance of this a priori enrichment of the reduced basis is illustrated on a typical aeronautical injector with Reynolds number of 45,000. This order reduction technique is able to recover large scale features for very long integration times (25 ms in the present case). Moreover, the reduced order model exhibits periodic fluctuations with a period of 2.2 ms corresponding to the time scale of the precessing vortex core (PVC) associated with this test case.

**An updated Gappy-POD to capture non-parameterized geometrical variation in fluid dynamics problems** In [ 4], a method is proposed to fill the gap within an incomplete turbulent and incompressible data field, while satisfying the topological and intensity changes of the fluid flow after a non-parameterized geometrical variation in the fluid domain. A single baseline large eddy simulation (LES) is assumed to be performed prior geometrical variations. The method enhances the Gappy-POD method proposed by Everson and Sirovich in 1995, in the case where the given set of empirical eigenfunctions is not sufficient and is not interpolant for the recovering of the modal coefficients for each Gappy snapshot by a least squares procedure. This is typically the case when we introduce non-parameterized geometrical modifications in the fluid domain. Here, after the baseline simulation, additional solutions of the incompressible Navier–Stokes equations are solely performed over a restricted fluid domain, that contains the geometrical modifications. These local LESs, called hybrid simulations, are performed by using the immersed boundary technique, which uses of a fluid boundary condition and the baseline velocity field. Then, the POD modes are updated using a local modification of the baseline POD modes in the restricted fluid domain. Furthermore, a physical correction of the latter enhanced Gappy-POD modal coefficients is obtained thanks to a Galerkin projection of the Navier–Stokes equations upon the new modes of the available data. This enhancement procedure on the global velocity reconstruction by the physical constraint was tested on a 3D semi-industrial test case of a typical aeronautical injection system and, a 2D laminar and unsteady incompressible test case. The speed-up relative to this new technique is equal to 100, which allows the fast exploration of two new designs of the aeronautical injection system.

### **6.2 Nonlinear Dimensionality Reduction via Auto-Encoder**

**Data-Targeted Prior Distribution for Variational AutoEncoder** In [ 1], Variational AutoEncoders (VAE) are used to study unsteady and compressible fluid flows in aircraft engines. Inferential methods enable the computation of sharp approximations of the posterior probability of the parameters of the VAE using the transient dynamics of the training velocity fields, and the generatation of plausible velocity fields. An important application is the initialization of such transient simulations. It is known by the Bayes theorem that the choice of the prior distribution is very important for the computation of the posterior probability, proportional to the product of likelihood with the prior probability. A new inference model is proposed, based on a new prior defined by the density estimate with the realizations of the kernel proper orthogonal decomposition coefficients of the available training data. It is illustrated that this inference model improves the results obtained with the usual standard normal prior distribution. This new generative approach can also be seen as an improvement of the kernel proper orthogonal decomposition method, for which we do not usually have a robust technique for expressing the pre-image in the input physical space of the stochastic reduced field in the feature high-dimensional space with a kernel inner product.

**A Bayesian Nonlinear Reduced Order Modeling Using Variational AutoEncoders** In [ 2], a new nonlinear projection based model reduction using convolutional VAEs. This framework is applied on transient incompressible flows. The accuracy is obtained thanks to the expression of the velocity and pressure fields in a nonlinear manifold maximising the likelihood on pre-computed data in the offline stage. A confidence interval is obtained for each time instant thanks to the definition of the reduced dynamic coefficients as independent random variables for which the posterior probability given the offline data is known. The parameters of the nonlinear manifold are optimized as the ones of the decoder layers of an autoencoder. The parameters of the conditional posterior probability of the reduced coefficients are the ones of the encoder layers of the same autoencoder. This reduced-order model is not a regression model over the offline pre-computed data. The numerical resolution of the ROM is based on the Chorin projection method. This new nonlinear projectionbased reduced order modeling is applied to a 2D Karman Vortex street flow and a 3D incompressible and unsteady flow in an aeronautical injection system.

**Deep multimodal autoencoder for crack criticality assessment** In continuum mechanics, the prediction of defect harmfulness requires to solve approximately partial differential equations with given boundary conditions. In [ 14], boundary conditions are learnt for tight local volumes (TLV) surrounding cracks in threedimensional volumes. A nonparametric data-driven approach is used to define the space of defects, by considering defects observed via X-Ray computed tomography. The dimension of the ambient space for the observed images of defects is huge. A nonlinear dimensionality reduction scheme is proposed in order to train a reduced latent space for both the morphology of defects and their local mechanical effects in the TLV. A multimodal autoencoder enables to mix morphological and mechanical data. It contains a single latent space, termed mechanical latent space. But this latent space is fed by two encoders. One is related to the images of defects and the other to mechanical fields in the TLV. The latent variables are input variables for a geometrical decoder and for a mechanical decoder. In this work, mechanical variables are displacement fields. The autoencoder on mechanical variables enables projectionbased model order reduction as proposed in the study of Lee and Carlberg. The main novelty of this work is a submodeling approach assisted by artificial intelligence. Here, for defect images in the test set, Dirichlet boundary conditions are applied to TLV. These boundary conditions are forecasted by the mechanical decoder with a latent vector predicted by the morphological encoder. For that purpose, a mapping is trained to convert morphological latent variables into mechanical latent variables, denoted "direct mapping". An "inverse mapping" is also trained for error estimation with respect to morphological predictions. Errors on mechanical predictions are close to 5% with simulation speed-up ranging for 3 to 120. It is illustrated that latent variables forecasted by the images of defects are prone to a better understanding of the predictions.

### **6.3 Piecewise Linear Dimensionality Reduction via Dictionary-Based ROM-Nets**

**Model order reduction assisted by deep neural networks (ROM-net)** In [ 12], a general framework is proposed for projection-based model order reduction assisted by deep neural networks. The proposed methodology, called ROM-net, consists in using deep learning techniques to adapt the reduced-order model to a stochastic input tensor whose nonparametrized variabilities strongly influence the quantities of interest for a given physics problem. In particular, the concept of dictionarybased ROM-nets is introduced, where deep neural networks recommend a suitable local reduced-order model from a dictionary. The dictionary of local reduced-order models is constructed from a clustering of simplified simulations enabling the identification of the subspaces in which the solutions evolve for different input tensors. The training examples are represented by points on a Grassmann manifold, on which distances are computed for clustering. This methodology is applied to an anisothermal elastoplastic problem in structural mechanics, where the damage field depends on a random temperature field. When using deep neural networks, the selection of the best reduced-order model for a given thermal loading is 60 times faster than when following the clustering procedure used in the training phase.

**Mechanical dissimilarity of defects in welded joints via Grassmann manifold and machine learning** Assessing the harmfulness of defects based on images is becoming more and more common in industry. At present, these defects can be

inserted in digital twins that aim to replicate in a mechanical model what is observed on a component so that an image-based diagnosis can be further conducted. However, the variety of defects, the complexity of their shape, and the computational complexity of finite element models related to their digital twin make this kind of diagnosis too slow for any practical application. In [ 18], a classification of observed defects enables the definition of a dictionary of digital twins. These digital twins prove to be representative of model-reduction purposes while preserving an acceptable accuracy for stress prediction. Nonsupervised machine learning is used for both the classification issue and the construction of reduced digital twins. The dictionary items are medoids found by a k-medoids clustering algorithm. Medoids are assumed to be well distributed in the training dataset according to a metric or a dissimilarity measurement. A new dissimilarity measurement between defects is proposed. It is theoretically founded according to approximation errors in hyper-reduced predictions. In doing so, defect classes are defined according to their mechanical effect and not directly according to their morphology. In practice, each defect in the training dataset is encoded as a point on a Grassmann manifold. This methodology is evaluated through a test set of observed defects totally different from the training dataset of defects used to compute the dictionary of digital twins. The most appropriate item in the dictionary for model reduction is selected according to an error indicator related to the hyper-reduced prediction of stresses. No plasticity effect is considered here (merely isotropic elastic materials), which is a strong assumption but which is not critical for the purpose of this work. In spite of the large variety of defects, accurate predictions of stresses for most of defects in the test set are obtained.

**Multimodal data augmentation for digital twining assisted by artificial intelligence in mechanics of materials** Digital twins in the mechanics of materials usually involve multimodal data in the sense that an instance of a mechanical component has both experimental and simulated data. These simulations aim not only to replicate experimental observations but also to extend the data. Whether spatially, temporally, or functionally, augmentation is needed for various possible uses of the components to improve the predictions of mechanical behavior. Related multimodal data are scarce, high-dimensional and a physics-based causality relation exists between observational and simulated data. In [ 5], a data augmentation scheme coupled with data pruning is proposed, in order to limit memory requirements for high-dimensional augmented data. This augmentation is desirable for digital twining assisted by artificial intelligence when performing nonlinear model reduction. Here, data augmentation aims at preserving similarities in terms of the validity domain of reduced digital twins. A specimen subjected to a mechanical test at high temperature is considered, where the as-manufactured geometry may impact the lifetime of the component. Hence, an instance is represented by a digital twin that includes 3D X-Ray tomography data of the specimen, the related finite element mesh, and the finite element predictions of thermo-mechanical variables at several time steps. There is, thus, for each specimen, geometrical and mechanical information. Multimodal data, which couple different representation modalities together, are hard to collect, and annotating them requires a significant effort. Thus, the analysis of multimodal data generally suffers from the problem of data scarcity. The proposed data augmentation scheme aims at training a recommending system that recognizes a category of data available in a training set that has already been fully analyzed by using high-fidelity models. Such a recommending system enables the use of a ROM-net for fast lifetime assessment via local reduced-order models.

**Real-Time Data Assimilation in Welding Operations Using Thermal Imaging and Accelerated Digital Twinning** Welding operations may be subjected to different types of defects when the process is not properly controlled and most defect detection is done a posteriori. The mechanical variables that are at the origin of these imperfections are often not observable in situ. In [ 16], an offline/online data assimilation approach is proposed, that allows for joint parameter and state estimations based on local probabilistic surrogate models and thermal imaging in real-time. Offline, the surrogate models are built from a high-fidelity thermomechanical finite element parametric study of the weld. The online estimations are obtained by conditioning the local models by the observed temperature and known operational parameters, thus fusing high-fidelity simulation data and experimental measurements.

**Mechanical assessment of defects in welded joints: morphological classification and data augmentation** In [ 15], a methodology is developed for classifying defects based on their morphology and induced mechanical response. The proposed approach is fairly general and relies on morphological operators and spherical harmonic decomposition as a way to characterize the geometry of the pores, and on the Grassman distance evaluated on FFT-based computations, for the predicted elastic response. The approach is implemented and detailed on a set of trapped gas pores observed in X-ray tomography of welded joints, that significantly alter the mechanical reliability of these materials. The space of morphological and mechanical responses is first partitioned into clusters using the "k-medoids" criterion and associated distance functions. Second, multiple-layer perceptron neuronal networks are used to associate a defect and corresponding morphological representation to its mechanical response. It is found that the method provides accurate mechanical predictions if the training data contains a sufficient number of defects representing each mechanical class. To do so, the original set of defects is supplemented by data augmentation techniques. Artificially-generated pore shapes are obtained using the spherical harmonic decomposition and a singular value decomposition performed on the pores signed distance transform.

**Data Augmentation and Feature Selection for Automatic Model Recommendation in Computational Physics** Classification algorithms have recently found applications in computational physics for the selection of numerical methods or models adapted to the environment and the state of the physical system. For such classification tasks, labeled training data come from numerical simulations and generally correspond to physical fields discretized on a mesh. Three challenging difficulties arise: the lack of training data, their high dimensionality, and the non-applicability of common data augmentation techniques to physics data. In [ 13], two algorithms are introduced to address these issues: one for dimensionality reduction via feature selection, and one for data augmentation. These algorithms are combined with a wide variety of classifiers for their evaluation. When combined with a stacking ensemble made of six multilayer perceptrons and a ridge logistic regression, they enable reaching an accuracy of 90% on the considered classification problem of nonlinear structural mechanics.

**Physics-informed cluster analysis and a priori efficiency criterion for the construction of local reduced-order bases** Nonlinear model order reduction has opened the door to parameter optimization and uncertainty quantification in complex physics problems governed by nonlinear equations. In particular, the computational cost of solving these equations can be reduced by means of local reduced-order bases. In [ 11], the benefits of a physics-informed cluster analysis for the construction of cluster-specific reduced-order bases are examined. The choice of the dissimilarity measure for clustering is fundamental and highly affects the performances of the local reduced-order bases. It is shown that clustering with an angle-based dissimilarity on simulation data efficiently decreases the intra-cluster Kolmogorov N-width. Additionally, an a priori efficiency criterion is introduced to assess the relevance of ROM-nets. This criterion also provides engineers with a very practical method for ROM-nets' hyperparameters calibration under constrained computational costs for the training phase. On five different physics problems, the physics-informed clustering strategy significantly outperforms classic strategies for the construction of local reduced-order bases in terms of projection errors.

### **6.4 Extension: Manifold Learning of Physics Problems Assisted by Black-Box Regressors**

**A priori compression of convolutional neural networks for wave simulators** Convolutional Neural Networks (CNNs) are seeing widespread use in a variety of fields, including image classification, facial and object recognition, medical imaging analysis, and many more. In addition, there are applications such as physics-informed simulators in which accurate forecasts in real time with a minimal lag are required. The current neural network designs include millions of parameters, which makes it difficult to install such complex models on devices that have limited memory. Compression techniques might be able to resolve these issues by decreasing the size of CNN models that are created by reducing the number of parameters that contribute to the complexity of the models. In [ 7], a compressed tensor format of convolutional layer is proposed a priori, before the training of the neural network, for finite element (FE) predictions of physical data. 3-way kernels or 2-way kernels in convolutional layers are replaced by one-way fiters. The overfitting phenomena will be reduced also. The time needed to make predictions or time required for training using the original CNN model would be cut significantly if there were fewer parameters to deal with. The priori compressed models is validated on physical data from a FE model solving a 2D wave equation. Tn the considered application, the proposed convolutional compression technique achieves equivalent performance in the prediction error as classical convolutional layers with fewer trainable parameters (around 20%) and lower memory footprint.

**Accelerated uncertainty quantification in impact simulations using generative adversarial networks and submodeling** The analysis of parametric and nonparametric uncertainties of very large dynamical systems requires the construction of a stochastic model of said system. Linear approaches relying on random matrix theory and principal component analysis can be used when systems undergo low-frequency vibrations. In the case of fast dynamics and wave propagation, a random generator of boundary conditions for fast submodels by using machine learning is investigated in [ 6]. It is illustrated that the use of non-linear techniques in machine learning and data-driven methods is highly relevant. Physics-Informed Neural Networks (PINNs) are a possible choice for a data-driven method to replace linear modal analysis. An architecture that supports a random component is necessary for the construction of the stochastic model of the physical system for non-parametric uncertainties, since the goal is to learn the underlying probabilistic distribution of uncertainty in the data. Generative Adversarial Networks (GANs) are suited for such applications, where the Wasserstein-GAN with gradient penalty variant offers improved convergence results for the considered problem. The objective of this approach is to train a GAN on data from a finite element method code so as to extract stochastic boundary conditions for faster finite element predictions on a submodel. The submodel and the training data have both the same geometrical support. It is a zone of interest for uncertainty quantification and relevant to engineering purposes. In the exploitation phase, the framework can be viewed as a randomized and parametrized simulation generator on the submodel, which can be used as a Monte Carlo estimator.

**MMGP: a Mesh Morphing Gaussian Process-based machine learning method for regression of physical problems under non-parameterized geometrical variability** When learning simulations for modeling physical phenomena in industrial designs, geometrical variabilities are of prime interest. While classical regression techniques prove effective for parameterized geometries, practical scenarios often involve the absence of shape parametrization during the inference stage, providing only mesh discretizations as available data. Learning simulations from such meshbased representations poses significant challenges, with recent advances relying heavily on deep graph neural networks to overcome the limitations of conventional machine learning approaches. Despite their promising results, graph neural networks exhibit certain drawbacks, including their dependency on extensive datasets and limitations in providing built-in predictive uncertainties or handling large meshes. In [10], a machine learning method that do not rely on graph neural networks is proposed. Complex geometrical shapes and variations with fixed topology are dealt with using well-known mesh morphing onto a common support, combined with classical dimensionality reduction techniques and Gaussian processes. The proposed methodology

can easily deal with large meshes without the need for explicit shape parameterization and provides crucial predictive uncertainties, which are essential for informed decision-making. In the considered numerical experiments, the proposed method is competitive with respect to existing graph neural networks, regarding training efficiency and accuracy of the predictions.

### **References**


**Open Access** This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.